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Objectives: 

The  objective  of  this  project  was  to  perform  a  series  of  systematic  numerical  experiments  in  order 
to  explore  new  concepts  for  active  control  of  turbulent  boundary  layers  for  viscous  drag  reduction. 
This  work  was  conducted  in  close  collaboration  with  AFOSRAJRI  Program  Integrated 
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Accomplishments: 

Three  different  control  schemes  have  been  developed.  The  first  control  scheme  utilized  an  adaptive 
controller  based  on  a  neural  network.  A  neural  network  was  constructed  and  trained,  and  then 
applied  to  turbulent  channel  flow  for  drag  reduction.  A  simple  control  network,  which  employs 
blowing  and  suction  at  the  wall  based  only  on  the  wall-shear  stresses  in  the  spanwise  direction,  was 
shown  to  reduce  the  skin  friction  by  as  much  as  20%  in  direct  numerical  simulations  of  a  low 
Reynolds  number  turbulent  channel  flow.  Also,  a  stable  pattern  was  observed  in  the  distnbution  of 
weights  associated  with  the  neural  network.  This  allowed  us  to  derive  a  simple  control  scheme 
that  produced  the  same  amount  of  drag  reduction.  This  simple  control  scheme  generates  optimum 
wall  blowing  and  suction  proportional  to  a  local  sum  of  the  wall-shear  stress  in  the  spanwise 
direction.  The  distribution  of  corresponding  weights  is  simple  ^d  localized,  thus  maWng  real 
implementation  relatively  easy.  The  result  of  this  work  was  published  in  Physics  of  Fluids  (Lee, 
Kim,  Babcock  and  Goodman,  Phys.  Fluids,  vol.  9,  no.  6, 1997). 

The  second  control  scheme  utilized  a  suboptimal  control  theory.  Two  simple  feedback  control  laws 
for  drag  reduction  were  derived  by  applying  a  suboptimal  control  theory  to  a  turbulent  chaniwl 
flow.  These  new  feedback  control  laws  required  pressure  or  shear-stress  information  only  at  the 
wall.  More  practical  control  laws  requiring  only  the  local  distribution  of  the  wall  pressure  or  one 
component  of  the  wall  shear  stress  were  also  derived  and  are  shown  to  work  equally  well.  The 
result  of  this  work  was  published  in  Journal  of  Fluid  Mechanics  (Lee,  Kim  and  Choi,  J.  Fluid 
Mech.,  vol.  338, 1998). 

The  third  one  deals  with  the  application  of  a  system  theory  approach  to  the  feedback  stabilization  to 
delay  the  transition  in  a  laminar  channel  flow.  It  was  shown  that  the  system  theory  approach  was  a 
valuable  tool  both  for  designing  control  systems  to  stabilize  the  flow  as  well  as  for  understanding 
the  physics  of  controlled  transitional  flows.  The  result  of  this  work  was  published  in  Journal  of 
Fluid  Mechanics  (Joshi,  Speyer  and  Kim,  J.  Fluid  Mech.,  vol.  332, 1997). 
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A  new  adaptive  controller  based  on  a  neural  network  was  constructed  and  applied  to  turbulent 
channel  flow  for  drag  reduction.  A  simple  control  network,  which  employs  blowing  and  suction  at 
the  wall  based  only  on  the  wall-shear  stresses  in  the  spanwise  direction,  was  shown  to  reduce  the 
skin  friction  by  as  much  as  20%  in  direct  numerical  simulations  of  a  low-Reynolds  number  turbulent 
channel  flow.  Also,  a  stable  pattern  was  observed  in  the  distribution  of  weights  associated  with  the 
neural  network.  This  allowed  us  to  derive  a  simple  control  scheme  that  produced  the  same  amount 
of  drag  reduction.  This  simple  control  scheme  generates  optimum  wall  blowing  and  suction 
proportional  to  a  local  sum  of  the  wall-shear  stress  in  the  spanwise  direction.  The  distribution  of 
corresponding  weights  is  simple  and  localized,  thus  making  real  implementation  relatively  easy. 
Turbulence  characteristics  and  relevant  practical  issues  are  also  discussed.  ©  1997  American 
Institute  of  Physics,  [81070-6631(97)02706-2] 


I.  INTRODUCTION 

The  ability  to  control  turbulent  flows  is  of  significant 
economic  interest.  Successful  control  of  turbulent  boundary 
layers  by  reducing  drag,  for  example,  can  result  in  a  substan¬ 
tial  reduction  in  operational  cost  for  commercial  aircraft  and 
marine  vehicles.  Recent  studies^’^  show  that  near-wall 
streamwise  vortices  are  responsible  for  high  skin-friction 
drag  in  turbulent  boundary  layers.  Some  attempts  have  been 
made  to  reduce  the  skin-friction  drag  by  controlling  the  in¬ 
teractions  between  these  vortices  and  the  wall.  Choi  et  al.^ 
for  example,  used  blowing  and  suction  at  the  wall  equal  and 
opposite  to  the  wall-normal  component  of  velocity  at 
y''"=  10.  They  showed  that  this  control  effectively  mitigated 
the  streamwise  vortices,  giving  approximately  25%  drag  re¬ 
duction  in  a  turbulent  channel  flow.  Although  the  method 
employed  in  their  work  is  impractical,  since  the  information 
at  y”''  =  10  is  usually  not  available,  it  demonstrates  a  control 
scheme  which  reduces  skin-friction  drag  by  manipulation  of 
the  near-wall  streamwise  vortices.  Another  way  to  control 
the  streamwise  vortices  is  to  impose  spanwise  oscillation  by 
a  moving  wall  or  externally  imposed  body  force.^’"^  These 
methods,  however,  require  a  large  amount  of  energy  input. 

A  systematic  approach  using  suboptimal  control  theory 
has  also  been  tried  in  the  past.  This  approach,  which  attempts 
to  minimize  a  cost  function,  was  successfully  applied  to  the 
stochastic  Burgers  equation.^  Moin  and  Bewley^  applied  a 
similar  approach  to  a  turbulent  channel  flow  to  achieve  up  to 
50%  drag  reduction.  The  control,  however,  requires  informa¬ 
tion  from  the  entire  velocity  field  inside  the  flow  domain  and 
excessive  computational  time,  making  it  impractical  to 
implement  in  real  situations.  For  practical  implementation,  a 
control  scheme  should  utilize  only  quantities  that  are  easily 

^Corresponding  author:  Telephone:  (310)825-4393;  Fax:  (310)206-4830; 

Electronic  mail:  jkim@turb.seas.ucla.edu 


measurable  at  the  wall,  and  should  be  fast  enough  to  be 
applied  in  real  time. 

The  objective  of  the  present  work  is  to  seek  wall  actua¬ 
tions,  in  the  form  of  blowing  and  suction  at  the  wall,  depen¬ 
dent  on  the  wall-shear  stress  to  achieve  a  substantial  skin- 
friction  reduction.  This  requires  knowledge  of  how  the  wall- 
shear  stresses  respond  to  wall  actuations,  i.e.,  the  correlation 
between  the  wall-shear  stresses  and  the  wall  actuations.  Be¬ 
cause  of  the  complexity  of  the  solutions  to  the  Navier- 
Stokes  equations,  however,  it  is  not  possible  to  find  such  a 
correlation  in  closed  form  or  to  approximate  it  in  simple 
form.  Instead,  we  use  a  neural  network  to  approximate  the 
correlation  which  then  predicts  the  optimal  wall  actuations  to 
achieve  the  minimum  value  of  the  skin-friction  drag.  Neural 
networks  have  been  used  to  obtain  complicated,  nonlinear 
correlations  without  a  priori  knowledge  of  the  system  that  is 
to  be  controlled.  Jacobson  and  Reynolds,^  for  instance,  used 
a  neural  network  to  obtain  about  7%  drag  reduction  in  their 
simulation  of  an  artificial  flow.  In  this  paper,  we  describe 
how  we  constructed  and  trained  a  neural  network  off-line, 
and  then  implemented  an  on-line  control  scheme  (blowing 
and  suction  at  the  wall)  for  drag  reduction  based  on  that 
neural  network.  Applying  this  control  scheme  to  direct  nu¬ 
merical  simulations  of  turbulent  channel  flow  at  low  Rey¬ 
nolds  number,  we  observed  about  20%  drag  reduction.  We 
then  describe  how  examination  of  the  weight  distribution 
from  the  on-line  neural  network  led  to  a  very  simple  control 
scheme  that  worked  equally  well  while  being  computation¬ 
ally  more  efficient. 

In  Sec.  II,  a  brief  description  of  the  architecture  of  the 
neural  network  used  in  the  present  work  is  given.  In  Sec.  HE, 
results  obtained  from  control  using  an  off-line  trained  net¬ 
work  are  presented,  while  results  obtained  from  control  using 
a  network  with  continuous  on-line  training  are  given  in  Sec. 
IV.  In  Sec.  V,  a  simple  control  law  based  on  the  weight 
distribution  from  the  successful  neural  network  control  is 
presented.  A  few  turbulence  statistics  are  given  in  Sec.  VI, 
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Sur^HCC  Shc^r  (dw/dy^) 
Stresses 

FIG.  1.  Neural  network  architecture. 


followed  in  Sec.  VII  by  a  discussion  of  practical  implemen¬ 
tation  and  the  conclusions. 

In  this  paper  we  use  (x,y,z)  for  the  streamwise,  wall- 
normal,  and  spanwise  coordinates,  respectively,  and 
(a,i;,w)  for  the  corresponding  velocity  components. 


II.  NEURAL  NETWORK 

In  this  section  we  describe  the  construction  of  a  neural 
network  to  learn  the  correlation  between  wall-shear  stresses 
and  wall  actuations  from  a  given  data  set.  Although  a  neural 
network  generally  requires  no  prior  knowledge  of  the  system 
(or  “plant”),  knowledge  about  the  near-wall  turbulence 
structures  provides  a  guideline  for  the  design  of  the  network 
architecture.  Initially  duldy  and  dwidy  at  the  wall  at  several 
instances  of  time  were  used  as  input  data  fields  and  the  ac¬ 
tuation  at  the  wall  was  used  for  the  output  data  of  the  net¬ 
work.  Experimentally  we  found  that  only  dwidy  at  the  wall 
from  the  current  instance  of  time  was  necessary  for  sufficient 
network  performance.  Because  we  wanted  the  output  to  be 
based  only  on  a  local  input  area,  we  designed  our  network 
using  shared  weights.  The  network  had  a  single  set  of 
weights  (a  template)  that  is  convolved  over  the  entire  input 
space  to  generate  output  values;  that  is,  we  used  the  same  set 
of  weights  for  each  data  point  and  the  training  involves  iter¬ 
ating  over  all  data  points.  The  template  extracts  spatially 
invariant  correlations  between  input  and  output  data.  The 
size  of  the  template  was  initially  chosen  to  include  informa¬ 
tion  about  a  single  streak  and  streamwise  vortex,  and  then 
was  varied  to  find  an  optimal  size. 

We  used  a  standard  two-layer  feedforward  network  with 
hyperbolic  tangent  hidden  units  and  a  linear  output  unit  (see 
Fig.  1).  The  functional  form  of  our  final  neural  network  is 


{N-\)n 


dw 


M+* 


(1) 


where  the  W’s  denote  weights,  N  is  the  total  number  of  input 
weights,  and  the  subscripts  j  and  k  denote  the  numerical  grid 
point  at  the  wall  in  the  streamwise  and  spanwise  directions 
respectively.  and  are  the  number  of  computational 


domain  grid  points  in  each  direction.  The  summation  is  done 
over  the  spanwise  direction.  Seven  neighboring  points 
(A— 7),  including  the  point  of  interest,  in  the  spanwise  di¬ 
rection  (corresponding  to  approximately  90  wall  units  with 
our  numerical  resolution)  were  found  to  provide  enough  in¬ 
formation  to  adequately  train  and  control  the  near-wall  struc¬ 
tures  responsible  for  the  high-skin  friction.  Note  that  the 
blowing  and  suction  are  applied  at  each  grid  location  accord¬ 
ing  to  (1)  as  a  numerical  approximation  of  distributed  blow¬ 
ing  and  suction  on  the  surface.  A  scaled  conjugate  gradient 
learning  algorithm^  was  used  to  produce  rapid  training.  For 
given  pairs  of  (vjl^ ,dw/dy\jk),  network  was  trained  to  mini¬ 
mize  the  sum  of  a  weighted-squared  error  given  by 

error=i2  E  (2) 

^  j  k 

where  is  the  desired  output  value  and  is  the  network 
output  value  given  by  Eq.  (1).  The  weights  were  initialized 
with  a  set  of  random  numbers.  Note  that  the  error  defined  in 
(2)  exponentially  emphasizes  (proportional  to  X)  large  actua¬ 
tions.  This  error  scaling  was  chosen  based  on  Choi  et  aV^ 
observation  that  large  actuations  are  more  important  for  drag 
reduction.  Usually  within  100  training  epochs,  the  error 
reached  its  asymptotic  limit. 

III.  OFF-LINE  TRAINING  AND  CONTROL 

As  an  initial  experiment  we  investigated  whether  we 
could  train  a  neural  network  to  predict  the  velocity  at 
>>■^  =  10  from  only  the  wall-shear  stresses.  The  rationale  be¬ 
hind  this  experiment  was  to  duplicate  an  existing  control  law 
using  only  measurements  available  at  the  wall  such  that  the 
output  from  the  network  could  be  used  as  input  to  the  actua¬ 
tor.  The  network  should  yield  a  similar  amount  of  drag  re¬ 
duction  to  that  obtained  by  Choi  et  al}  through  prediction  of 
the  velocity  at  3^*^  =  10.  The  training  data  consisted  of  100 
time  steps  of  output  obtained  from  a  numerical  simulation  of 
channel  flow  under  the  control  employed  by  Choi  et  al^  i.e., 
using  the  wall-normal  velocity  at  y '*'  =  10.  The  flow  regime  is 
turbulent  channel  flow  with  Re^=  100,  where  Re^  is  the  Rey¬ 
nolds  number  based  on  the  wall-shear  velocity,  m^,  and  the 
channel  half-width,  5.  All  numerical  simulations  presented  in 
this  paper  were  obtained  using  a  modified  version  of  Kim 
etaV^  spectral  code  with  the  computational  domain 
(47r,2,47r/3)5,  and  a  grid  resolution  of  (32,  65,  32)  in  the 
(x,y,z)  directions,  respectively.  Each  time  step  contained  a 
32X32  array  of  input  values  {dwldy\y^)  and  corresponding 
actuations,  —v  aty‘*'  =  10. 

We  trained  several  networks  with  one  hidden  unit  and 
different  sized  input  templates:  7X1,  7X3,  7X5,  9X1, 
9X3,  llXl,  and  11X3  (the  number  of  input  points  in  the 
spanwise  direction  by  the  number  of  input  points  in  the 
streamwise  direction).  After  training  was  completed,  the  dis¬ 
tribution  of  input  weights  was  examined  to  see  whether  there 
was  a  discernible  pattern  in  the  input  template.  The  weight 
distributions  in  the  spanwise  direction  at  the  same  stream- 
wise  location  for  different  input  templates  are  shown  in  Fig. 
2.  The  same  pattern  for  all  7  input  template  sizes  was  ob¬ 
served,  and  is  similar  to  a  finite  difference  approximation  of 
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FIG.  2.  Weight  distribution  from  off-line  trainings  for  various  sizes  of  the 
input  template:  7  X  1,7  X  3,7  X  5,9  X  1,9  X  3, 11  X  1, 11  X  3.Weightsare 
normalized  by  Wp 


spanwise  differentiation.  Increasing  the  number  of  hidden 
units  only  marginally  improved  performance,  whereas  in¬ 
creasing  the  template  size  significantly  reduced  the  final 
training  error. 

We  then  applied  a  control  scheme  to  a  regular  channel 
flow  by  fixing  the  final  input  weights  obtained  from  the  off¬ 
line  training.  This  was  perhaps  a  somewhat  naive  approach 
since  the  weights  were  obtained  from  fully  controlled  flow 
data  and  the  shear  stresses  used  for  the  training  were  already 
altered  by  the  actuations.  Nevertheless,  two  cases  were 
tested:  a  control  scheme  based  on  7  weights  in  the  spanwise 
direction,  and  another  based  on  the  same  7  weights  plus  3 
more  in  the  immediate  downstream  location.  These  10 
weights  were  chosen  because  among  all  weights  obtained 
from  the  off-line  training  they  had  non-negligible  values.  In 
Fig.  3,  the  mean  shear  stress  variations  at  the  wall  (i.e.,  drag) 
obtained  with  these  two  controls  are  plotted  along  with  the 
no-control  case.  Nearly  18%  drag  reduction  was  achieved, 
with  slightly  better  performance  from  the  10-weight  net¬ 
work. 

These  results  demonstrate  that  a  correlation  exists  be¬ 
tween  the  shear  stresses  at  the  wall  and  the  desired  actua- 


FIG.  3.  Mean  wall-shear  stress  histories  for  various  control  laws  compared 
to  the  no-control  case:  — ,  no  control;  — control  with  7  fixed  weights 
obtained  from  off-line  control;  control  with  10  fixed  weights  from  the 
same  off-line  control.  All  runs  with  control  shown  in  this  paper  were  inte¬ 
grated  over  at  least  20  time  units  that  is  much  longer  than  the  typical  turn¬ 
over  time  of  the  streamwise  vortices,  1  time  unit. 
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FIG.  4.  Schematic  representation  of  adaptive  inverse  model  control. 


tions,  and  that  control  based  on  this  correlation  produces  a 
significant  amount  of  drag  reduction.  This  fixed-weight  con¬ 
trol  scheme,  however,  was  deduced  from  fully  controlled 
flow  data;  thus  it  does  not  guarantee  the  same  performance 
for  other  flow  situations. 


IV.  ON-LINE  CONTftOL 


In  the  previous  section  we  showed  how  successful  con¬ 
trol  based  on  off-line  training  can  be  obtained.  However, 
since  the  system  we  are  trying  to  control  is  time- varying  and 
nonlinear,  this  approach  is  not  likely  to  generalize  well.  Con¬ 
tinuous  on-line  training  allows  a  controller  to  adapt  to  the 
evolution  of  the  system.  In  this  section,  we  describe  an  adap¬ 
tive  controller  for  on-line  training  and  control. 

There  are  various  schemes  for  on-line  neural  network 
control.  The  most  direct  scheme  is  adaptive-inverse  model 
control.  A  schematic  representation  of  this  approach  is 
shown  in  Fig.  4.  Here  the  plant  denotes  the  numerical  solver 
of  the  Navier-Stokes  equations.  This  configuration  employs 
a  neural  network  to  model  the  (possibly  time-varying)  in¬ 
verse  plant  mapping  from  wall-shear  stress  dwldy\y^  to  the 
wall-normal  actuations,  and  then  uses  a  copy  of  the  model  as 
the  controller  with  the  desired  shear  stresses  as  input.  One 
restriction  of  this  technique  is  that  it  usually  requires  an  ini¬ 
tial  model  training  phase  using  random  plant  inputs  and  cor¬ 
responding  plant  outputs.  This,  however,  caused  no  serious 
problems  since  usually  one  time  step  was  enough  for  model 
training  because  of  the  shared- weight  architecture  in  our  ap¬ 
plication  (each  time  is  1024  data  points  and  our  networks 
only  have  a  few  weights).  Once  the  model  represents  a  rea¬ 
sonably  close  approximation  to  the  actual  plant  inverse,  a 
copy  is  then  implemented  as  a  feedforward  controller. 

The  desired  inputs  to  the  controller  are  a  fractional  re¬ 
duction  in  the  shear  stress  from  the  previous  step,  i.e.. 


(3) 


where  0<  7)<  1 .  Based  on  the  networks  and  techniques  we 
investigated,  this  indirect  suppression  of  dwldy\y^ ,  instead  of 
duldy\^,  turns  out  to  be  more  efficient  in  achieving  drag 
reduction.  The  output  of  the  controller,  which  is  the  input  to 
the  plant,  is  the  predicted  actuation  necessary  to  produce  this 
shear-stress  reduction.  The  quantity  7]  should  be  chosen  such 
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FIG.  5.  Distribution  of  the  weights:  0,  from  on-line  training;  — , 
A(1  -cos(-n7))//  ^2"^=  13. 


that  sets  of  the  desired  large  amplitude  outputs  are  among  the 
sets  that  are  well  represented  by  the  training  sets.  Good  per¬ 
formance  was  achieved  for  the  range  of  77=0.8-0.85. 

A  turbulent  channel  flow  at  Reynolds  number 
Rey=  100  was  used  to  test  the  neural  network.  We  allowed 
all  the  weights  in  the  network  to  adapt  and  examined  the 
input  template  pattern  after  each  time  step.  As  the  control 
began,  the  weight  distribution  immediately  assumed  a  fixed 
pattern  similar  to  the  off-line  pattern.  There  was  no  appre¬ 
ciable  change  in  the  relative  magnitudes  of  the  template 
weights  over  time,  indicating  that  the  pattern  is  preserved. 
The  absolute  magnitudes,  however,  did  vary  indicating  the 
need  for  gain  and  bias  adaptation  for  each  layer.  The  number 
of  hidden  layers,  the  number  of  the  hidden  units,  the  size  of 
the  input  template,  and  the  error  scale  of  the  training  error 
[X  of  Eq.  (2)]  were  1,  1,  7X1,  and  5,  respectively.  We 
varied  the  values  of  the  error  scale  and  found  5  to  be  opti¬ 
mum,  with  larger  values  causing  an  instability  in  training. 
We  also  varied  the  number  of  the  hidden  units,  but  the  net¬ 
work  with  a  single  hidden  unit  produced  the  best  result.  The 
converged  weight  distributions  for  20  consecutive  time  steps 
after  an  asymptotic  state  was  reached  are  shown  in  Fig.  5. 
The  pattern  is  only  slightly  different  from  the  one  obtained 
from  the  off-line  training  (see  Fig.  2).  It  should  be  noted  that 
this  pattern  emerged  immediately  after  the  on-line  control 
began. 

Time  histories  of  the  wall-shear  stress  for  3  different 
input  template  sizes  are  shown  in  Fig.  6.  All  other  network 
parameters  are  kept  the  same.  The  abscissa  represents  non- 
dimensional  time  normalized  by  and  S.  One  time  unit 
roughly  corresponds  to  the  time  for  a  fluid  particle  traveling 
with  the  centerline  velocity  to  move  about  20  channel  half 
widths.  The  computations  were  carried  out  long  enough  to 
ensure  that  a  statistically  steady  state  was  reached.  As  the 
control  began,  the  drag  quickly  drops  to  about  80%  of  that 
observed  without  control,  for  the  two  cases  with  template 
sizes  of  7  X 1  and  9X1.  The  template  size  of  5  X 1 ,  however, 
did  not  produce  as  much  reduction,  indicating  that  at  least  7 
spanwise  points  (which  extended  about  90  wall  units,  with 
our  grid  resolution)  should  be  used  for  good  performance.  At 
the  initial  stage  of  our  study,  we  tested  a  control  using  both 
^w/^y\y^^  and  du/dy\^  as  input,  with  a  7X5  input  template. 
It  produced  about  the  same  amount  of  reduction,  but  with 
excessive  training  time  due  to  a  significantly  larger  number 


FIG.  6.  Mean  wall-shear  stress  histories  for  on-line  control  with  different 
input  template  sizes:  — ,  no  control;  control  with  5X1;  — ,  control  with 
7X1;  — control  with  9X 1, 


of  weights.  Furthermore,  it  did  not  produce  coherent  patterns 
in  the  weight  distributions. 

Since  the  7  weights  showed  a  fixed  pattern,  we  fixed  the 
input  template  weights  to  this  pattern  and  used  a  single  hid¬ 
den  unit  network,  giving  only  4  adaptable  parameters  (a  bias 
and  gain  for  each  layer).  Iliis  simplified  network  had  the 
following  functional  form: 

Vjk  =  tanh(  Wi,g  -  W^)  -  (4) 

with 


^  Sw 


j.k+i 


(5) 


where  W/’s  are  the  fixed-weight  pattern  obtained  from  the 
previous  on-line  control.  On-line  control  using  this  network 
produced  a  similar  amount  of  drag  reduction.  The  weight 
variations  with  time  were  also  monitored.  The  bias  weights 
(Wc  and  W^)  were  negligibly  small.  Those  controlling  the 
gain  (Wa  and  W*),  which  had  finite  values,  changed  in  time 
significantly,  although  the  product  of  the  two  gain-weights 
remained  almost  constant.  This  suggests  that  effective  con¬ 
trol  can  be  achieved  by  simply  using  g  with  an  adjustable 
amplitude.  This  will  be  discussed  in  the  following  section. 


V.  A  SIMPLE  CONTROL  SCHEME 


Since  the  control  based  on  the  network  given  by  Eq.  (4) 
produced  a  substantial  reduction  in  drag,  this  section  devel¬ 
ops  a  control  scheme  based  only  on  the  weighted  sum  of  the 
wall-shear  stress.  The  distribution  of  the  weights  can  be  ap¬ 
proximated  by  (see  Fig.  5): 


l-cos(^;) 

W;  =  A - ^ - 

^  J 


(6) 


where  y  =  0  corresponds  to  the  point  where  the  control  is 
applied.  Since  only  the  relative  values  are  important,  the  con¬ 
stant  A  has  no  special  meaning.  A  computation  with  a  twice 
higher  resolution  was  run  to  confirm  that  Eq.  (6)  gives  the 
proper  form  of  the  weight  distribution  (see  Fig.  7).  It  turns 
out  that  Eq.  (6)  is  the  inverse  Fourier  transform  of  ik^/\k^\y 
where  is  the  wave  number  in  the  spanwise  direction,  for 
the  finite  maximum  wave  number,  i.e.. 
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FIG.  7.  Distribution  of  the  weights:  0,  from  on-line  training; — , 
A  ( 1  —  cos(TTj))/j.  Az  ■*■  =  6 . 


FIG.  9.  Mean  wall-shear  stress  histories  for  various  control  schemes  com¬ 
pared  to  the  no-control  case:  — ,  no  control;  — ,  on-line  control  with  neural 
network  with  7X1  template;  •••  control  with  7  fixed  weights. 


jj^cxp{-ikzZ)dk^^2 


l-cos{k^z) 


(7) 


where  kn—'rrfAz  is  the  maximum  wave  number  and  Az  is 
the  numerical  grid  spacing.  Replacing  z  with  jAz  in  the 
right-hand~side  of  Eq.  (7)  leads  to  Eq.  (6).  From  this  result 
and  the  convolution  theorem,  one  can  suggest  the  following 
simple  control  law  for  wall  transpiration: 


ik^  dw  I  d  I  dw\ 

^  ^z\dyj 


(8) 


where  the  “hat”  denotes  a  Fourier  transformed  quantity  and 
C  is  a  positive  scale  factor  determining  the  amplitude  of  the 
actuation.  Equation  (8),  which  should  produce  similar  drag 
reduction  as  the  result  with  7  fixed  weights  [Eq.  (4)],  implies 
that  the  optimum  blowing  and  suction  at  the  wall  is  propor¬ 
tional  to  d{dwl dy)l dz  with  the  high  wave  number  compo¬ 
nent  suppressed  by  y\k^>  Note  that  the  wall-normal  actua¬ 
tion  proportional  to  ^(^w/^y)/^z  counteracts  the  up-and- 
down  motion  induced  by  a  streamwise  vortex  (see  Fig.  8), 
consistent  with  that  used  by  Choi  etal}  Smith  etal}^  dis¬ 
cussed  the  stability  of  the  near- wall  turbulence  structure  to  a 
local  adverse  pressure  gradient  imposed  near  the  surface. 


FIG.  8.  A  schematic  of  a  flow  field  induced  by  a  streamwise  vortex  and 
corresponding  d{dwl dy\^)! dz. 


The  suction  at  the  left  side  of  the  vortex  in  Fig.  8  can  relieve 
such  adverse  pressure  gradient,  thus  preventing  the  surface- 
layer  eruption.  Choi  et  al}  reported  that  actuation  based  only 
on  the  local  value  of  d{dwldy)ldzy  which  is  also  a  dominant 
term  of  the  Taylor  expansion  of  v  at  some  distance  away 
fi*om  the  wall,  did  not  produce  such  a  substantial  drag  reduc¬ 
tion.  This  indicates  that  suppression  of  high  wave  number 
components  by  or  equivalently  using  locally  weighted 
value  of  dwidy,  plays  a  key  role  in  reducing  drag.  Note  that 
the  weights  at  even  numbered  grid  points  away  from  the 
center  point  vanish.  This  is  beneficial  for  physical  implemen¬ 
tation,  since  sensors  and  actuators  cannot  be  placed  at  the 
same  location.^^  Because  the  Fourier  integral  is  computed  for 
a  finite  value  of  k^^ ,  the  values  of  the  weight  at  noninteger 
j  in  Eq.  (6)  have  no  meaning.  Equation  (8)  is  equivalent  to 


(Ar-i)/2  . 

Vjk=C  s  W,— 


7,it+i 


(9) 


where  Wi  is  given  by  Eq.  (6).  The  magnitude  of  the  weights 
decays  with  increasing  distance  from  the  center,  which  al¬ 
lows  for  good  approximation  using  only  a  small  number  of 
weights,  i.e.,  successful  control  requires  only  local  values  of 
the  shear  stress.  A  natural  concern  is  how  different  grid  spac- 
ings  affect  the  control.  Since  the  above  control  law  [Eq.  (9)] 
is  simply  another  expression  of  Eq.  (8),  which  is  a  good 
approximation  as  long  as  k^^  is  large  enough,  the  control 
should  be  relatively  independent  of  resolution.  This  was  con¬ 
firmed  by  a  computation  with  a  higher  resolution. 

Control  based  on  Eq.  (9)  with  7  points  produced  the 
same  amount  of  drag  reduction  (20%)  as  the  neural  network 
control  (see  Fig.  9).  The  constant  C  is  chosen  so  that  the 
root-mean-squared  (rms)  value  of  the  actuation  is  kept  at 
O.ISm^.  Blowing  and  suction  of  this  magnitude  at  the  wall 
suppress  the  near-wall  streamwise  vortices  by  counteracting 
the  up-and-down  motions  associated  with  these  vortices  (see 
Fig.  8). 

The  fact  that  a  control  scheme  that  is  essentially  linear 
[Eq.  (9)]  produced  the  same  amount  of  drag  reduction  as 
on-line  control  with  a  nonlinear  network  led  to  the  following 
question:  Can  the  same  weight  pattern  be  captured  when 
on-line  control  with  a  “linear”  network  is  used?  To  answer 
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FIG.  10.  Contours  of  streamwise  vorticity  in  a  cross-flow  plane:  (a)  no 
control;  (b)  control  using  7  fixed  weights.  The  contour  level  increment  is  the 
same  for  both  figures.  Negative  contours  are  chain-dotted. 


this  question,  we  tested  an  on-line  controller  utilizing  a  linear 
network.  The  hyperbolic  tangent  function  in  Eq.  (1)  was  re¬ 
placed  with  a  linear  function  and  the  hidden  layer  was  re¬ 
moved.  All  other  parameters  were  kept  the  same  as  the  non¬ 
linear  network  case.  The  linear  network  seemed  to  produce 
almost  the  same  pattern  in  weights  as  the  one  produced  by 
our  nonlinear  network.  Also,  drag  was  reduced  by  the  same 
amount  The  pattern  that  was  obtained  by  averaging  over 
time,  however,  was  not  well  preserved  in  time.  The  weight 
of  Eq.  (1),  for  example,  frequently  showed  order  of  mag¬ 
nitude  excursions  resulting  in  a  large  standard  deviation,  1.5, 
compared  to  the  value  of  the  weight,  0.3  (weights  were  nor¬ 
malized  by  the  maximum  amplitude  Wi),  The  nonlinear  net¬ 
work,  however,  produced  more  constant  weights;  the  stan¬ 
dard  deviation  of  was  only  0.15.  This  result  indicates  that 
the  nonlinearity  should  add  robustness  by  bounding  the 
weights.  This  can  simplify  hardware  implementation  by  lim¬ 
iting  the  signals  and  subsequently  the  necessary  dynamic 
range  of  the  processing  circuitry. 


VI.  TURBULENCE  STATISTICS 

The  computed  flow  fields  for  a  no-control  case  and  a 
successful  control  case,  based  on  the  7-point  weighted  sum 
of  dw/dy\^  [Eq.  (9)],  were  examined  to  investigate  the 
mechanism  by  which  the  drag  reduction  is  achieved.  The 
most  salient  feature  of  the  controlled  case  was  that  the 
strength  of  the  near-wall  streamwise  vortices  was  drastically 
reduced.  In  Fig.  10,  contours  of  streamwise  vorticity  in  a 
cross  plane  are  shown.  This  result  further  substantiates  the 
notion  that  a  successful  suppression  of  the  near-wall  stream- 


FIG.  11.  Probability-density  function  of  the  wall-shear  stress:  — .  no  con¬ 
trol;  _ ,  control  with  7  fixed  weights.  The  area  under  each  curve  is  normal¬ 

ized  to  one. 

wise  vortices  leads  to  a  significant  reduction  in  drag.  Note 
that  for  the  controlled  case  the  wall  actuations  were  applied 
at  both  walls. 

The  probability-density  function  of  the  wall-shear  stress 
in  the  streamwise  direction  is  shown  in  Fig.  11.  It  is  evident 
that  the  control  is  very  effective  in  suppressing  large  fluctua¬ 
tions,  thus  reducing  the  mean-skin  friction.  Furthermore,  the 
root-mean-square  (rms)  values  of  turbulent  fluctuations  in 
the  wall  region  are  also  reduced  as  shown  in  Fig.  12(a).  The 
same  trend  was  observed  by  Choi  et  dl?"  The  rms  vorticity 
fluctuations  are  also  significantly  reduced,  except  for  very 


FIG.  12.  Root-mean-square  fluctuations  normalized  by  the  wall  v^abfe 
— ,  no  control;  control  with  7  fixed  weights,  (a)  Velocity  fluctuations;  (b) 
vorticity  fluctuations. 
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FIG.  13.  Contours  of  the  wall  actuation:  (a)  control  using  dy^fdy\y,  with  7 
fixed  weights;  (b)  control  using  information  at  y"**- 10.  Negative  contours 
are  chain-dotted. 


close  to  the  wall.  The  increase  for  the  latter  is  caused  by 
additional  dvidz  due  to  the  wall  actuation.  The  rms  fluctua¬ 
tions  of  at  the  wall,  which  is  mainly  due  to  duldy  at  the 
wall,  is  also  decreased.  This  indicates  that  the  control  scheme 
led  to  a  reduction  in  the  mean  shear  and  its  variance  at  the 
wall  by  suppressing  large  fluctuations  as  shown  in  Fig.  11. 
The  reduction  in  the  rms  fluctuations  in  Wj,.  and  cOy  indicates 
that  the  control  scheme  indeed  reduced  Ae  strength  of  the 
near-wall  streamwise  vortices  and  the  wall-layer  streaks. 

Figure  13  compares  the  distribution  of  wall  actuation 
used  in  our  control  with  that  from  Choi  etaV^  u -control 
using  the  information  at  y'’‘=10  for  the  same  wall-shear 
stress  distribution.  The  corresponding  wall-shear  stress  dis¬ 
tribution  is  also  shown  in  Fig.  13.  The  wall  actuations  indi¬ 
cate  a  strikingly  similar  distribution  to  each  other,  even 
though  the  wall  actuation  of  our  control  is  based  only  on  the 
wall-shear  stress  dwldy\y, .  Basically  our  control  scheme  de¬ 
tects  edges  of  locally  high-shear  stress  regions  by  measuring 
the  spanwise  variation  of  dwidy,  and  applies  appropriate 
wall  actuation  as  shown  in  Figs.  13(a)  and  13(b).  Since  high- 
shear  stress  regions  are  usually  elongated  in  the  streamwise 
direction,  only  spanwise  variation  is  necessary  for  detecting 
the  edges.  Since  dwidy  is  a  direct  measurement  of  stream- 
wise  vortices  (see  Fig.  8),  it  provides  more  appropriate  infor¬ 
mation  than  duldy.  This  is  consistent  with  our  finding  that 
dwfdy  at  several  points  in  spanwise  direction  are  enough  for 
good  performance  of  control. 

VII.  DISCUSSION  AND  CONCLUSION 

We  have  presented  a  successful  application  of  a  neural 
network  to  turbulence  control  for  drag  reduction.  First  we 


were  able  to  construct  and  train  a  neural  network  off-line  to 
find  a  correlation  between  the  wall-shear  stress  and  the  de¬ 
sired  wall  actuations  coming  from  the  information  at 
y  +  =  10.  Based  on  the  optimal  network  structure  from  the 
off-line  training,  we  successfully  implemented  an  on-line  in¬ 
verse  model  controller  in  numerical  experiments  of  a  turbu¬ 
lent  channel  flow,  resulting  in  about  20%  drag  reduction. 
Finally,  we  were  able  to  deduce  a  simple  control  scheme 
[Eq.  (9)]  for  drag  reduction  based  on  observation  of  the 
weight  distribution  from  a  successful  control  case.  This  con¬ 
trol  scheme  uses  a  minimal  amount  of  wall-shear  stress  in¬ 
formation  and  requires  only  simple  operations,  thus  render¬ 
ing  a  scheme  whose  actual  hardware  implementation  would 
be  relatively  easy. 

Beginning  with  a  nonlinear  network,  we  found  a  simple 
linear  control  scheme  from  the  pattern  of  the  weights,  sug¬ 
gesting  the  use  of  a  linear  network.  Although  a  linear  net¬ 
work  produced  a  similar  pattern,  the  pattern  was  not  well 
preserved  in  time.  The  weights  can  vary  significantly  in¬ 
creasing  the  difficulty  in  circuit  implementation  due  to  the 
larger  required  dynamic  range.  The  nonlinear  network,  how¬ 
ever,  produced  bounded  weights  simplifying  hardware 
implementation.  This  suggests  that  a  certain  amount  of  non¬ 
linearity  is  beneficial  to  capture  a  stable  coherent  pattern  for 
our  system. 

Although  the  present  control  scheme  is  a  significant  im¬ 
provement  over  earlier  approaches  that  require  velocity  in¬ 
formation  within  the  flow,  there  is  a  technical  issue  before 
the  present  control  scheme  can  be  implemented  in  real  prac¬ 
tice.  Precise  control  of  blowing  and  suction  distributed  over 
a  surface  is  difficult  to  implement.  Other  approaches,  such  as 
surface  movements  by  deformable  surface,  may  prove  to  be 
more  practical. 

There  are  also  several  fundamental  issues  that  must  be 
addressed.  First,  our  numerical  experiments  were  performed 
in  a  very  low  Reynolds  number  flow.  It  remains  to  be  seen 
whether  the  same  control  scheme  extends  to  higher  Reynolds 
numbers.  If  the  main  cause  of  high-skin  friction  in  turbulent 
boundary  layers  at  higher  Reynolds  numbers  is  also  due  to 
the  near-wall  streamwise  vortices,  the  same  scheme  should 
work  equally  well.  Detection  of  these  vortices  through  the 
wall-shear  stress  would  become  more  difficult,  however,  be¬ 
cause  the  scales  associated  with  these  vortices  decrease  as 
the  Reynolds  number  increases. 

Another  important  issue  is  the  influence  of  the  spatial 
resolution  of  sensors  and  actuators  on  performance.  We 
showed  that  a  simple  control  scheme  with  4  neighboring 
points  in  the  spanwise  direction  (3  among  7  weights  are  ze¬ 
ros),  which  span  90  wall  units,  performed  very  well.  This 
suggests  that  the  distance  between  sensors  should  be  about 
20  wall  units.  As  the  Reynolds  number  increases,  the  physi¬ 
cal  separation  between  sensors  must  decrease.  We  note  in 
passing  that  recent  advances  in  micromachined  sensors  and 
actuators  make  these  scales  feasible,'^  and  the  present  work 
is  part  of  a  joint  research  project  aimed  at  integrating  micron¬ 
sized  sensors  and  actuators  with  analog  VLSI  control  cir¬ 
cuitry. 

A  third  issue  is  determining  an  appropriate  scale  factor 
C  in  Eq.  (9).  We  tested  a  range  of  C  and  empirically  found 
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that  values  yielding  an  actuation  rms  value  between  O.lw,. 
and  0.15m^  gave  the  best  performance.  Smaller  values  re¬ 
sulted  in  less  reduction,  while  larger  ones  caused  rapid  fluc¬ 
tuations  of  the  wall-shear  stress  over  time.  The  experiments 
were  only  performed  for  Re^=l(X)  and  180,  from  which  the 
optimum  value  was  deduced,  but  we  expect  that  the  same 
amplitude  level  should  produce  similar  reductions  for  higher 
Reynolds  number  flows.  We  further  note  that  for  this  ampli¬ 
tude  the  required  power  input  per  unit  area  to  produce  the 
actuation,  p„v„+0.5pvl,(~0.lpul)  was  significantly  less 
compared  to  the  power  saved  due  to  the  reduced  drag, 
&,CffCfr^Uc('-3,5pul)  for  20%  reduction  in  Cf,  Here 
Pw  >  P>  '^w  and  Uc  are  the  wall  pressure,  density,  friction 
coefficient,  averaged  wall-shear  stress  in  the  streamwise  di¬ 
rection,  and  centerline  velocity,  respectively.  Note  that  we 
did  not  account  for  device  loss  to  deliver  blowing  and  suc¬ 
tion  nor  the  power  consumed  by  sensors  in  this  estimate. 

One  final  practical  issue  worth  mentioning  is  the  time 
delay  between  sensing  and  actuation.  None  was  included  in 
any  of  our  numerical  experiments.  In  a  real  situation,  how¬ 
ever,  there  will  be  a  finite  time  delay  between  sensing  and 
actuation.  The  process  time  for  the  simple  control  scheme 
should  be  negligible  since  a  weighted  summation  is  an  in¬ 
volved  process,  thus  not  imposing  any  restriction  on  the  time 
delay.  Ideally,  sensors  should  register  the  response  of  the 
near-wall  structures  to  be  controlled  due  to  wall  actuations, 
which  will  take  a  certain  amount  of  time.  Blowing  and  suc¬ 
tion  at  the  wall,  however,  can  produce  an  immediate  spurious 
response  at  nearby  sensors  due  to  local  conservation.  One 
possible  remedy  to  this  problem  is  an  under-relaxation  of  the 
actuation  signal  with  past  signals,  accounting  for  the  re¬ 
sponse  of  a  longer  time  scale.  For  example,  an  underrelax¬ 
ation  scheme  using  the  following  formula 


^jk 


j*k  +  i 


(10) 


could  take  into  account  all  past  signals  with  a  single  param¬ 
eter  f  (0<^<  1).  We  used  this  approach  successfully  in  our 
numerical  experiment  for  Re,^"  found  that  there  is 

an  optimum  ^  depending  on  temporal  resolution. 

The  most  significant  finding  of  the  present  study  is  that  a 
single  span  wise  strip  of  is  enough  to  achieve  sig¬ 

nificant  drag  reduction.  Additional  information  about 
^w/^y\y^;  in  the  streamwise  direction  or  du/dy\^  in  either 
direction  reduced  the  efficiency  of  our  neural-network  based 
control.  We  also  tested  different-sized  input  templates  and 
found  that  as  the  template  size  increased  in  the  streamwise 
direction,  the  training  time  became  excessive  and  the  fixed 
pattern  in  the  weight  disappeared,  indicating  that  too  many 
unnecessary  weights  caused  training  problems. 
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Two  simple  feedback  control  laws  for  drag  reduction  are  derived  by  applying  a 
suboptimal  control  theory  to  a  turbulent  channel  flow.  These  new  feedback  control 
laws  require  pressure  or  shear-stress  information  only  at  the  wall,  and  when  applied 
to  a  turbulent  channel  flow  at  Re.^  =  110,  they  result  in  16-22%  reduction  in  the 
skin-friction  drag.  More  practical  control  laws  requiring  only  the  local  distribution  of 
the  wall  pressure  or  one  component  of  the  wall  shear  stress  are  also  derived  and  are 
shown  to  work  equally  well. 


1.  Introduction 

Recent  studies  have  shown  that  near-wall  streamwise  vortices  are  responsible 
for  high  skin-friction  drag  in  turbulent  boundary  layers.  Many  attempts  aiming  at 
controlling  these  vortices  have  been  made  in  order  to  achieve  a  skin-friction  drag 
reduction  in  turbulent  boundary  layers.  Most  of  such  attempts,  however,  have  been 
ad  hoc,  largely  based  on  physical  intuition.  The  active  control  by  Choi,  Moin  &  Kim 
(1994),  for  example,  used  blowing  and  suction  at  the  wall  that  is  equal  and  opposite 
to  the  wall-normal  component  of  the  velocity  at  y"*"  =  10,  resulting  in  as  much  as 
25%  reduction  in  their  numerical  simulation.  This  approach,  however,  is  impractical 
since  the  required  velocity  information  at  y'*'  =  10  is  not  normally  available.  For 
any  practical  implementation  a  control  scheme  should  be  based  solely  on  quantities 
measurable  at  the  wall. 

A  more  systematic  approach  based  on  an  optimal  control  theory  was  proposed 
by  Abergel  &  Temam  (1990).  Choi  et  al.  (1993)  proposed  a  ‘suboptimal’  control 
procedure,  in  which  the  iterations  required  for  a  global  optimal  control  were  avoided 
by  seeking  an  optimal  condition  over  a  short  time  period.  The  suboptimal  control 
procedure  was  successfully  applied  to  control  of  the  Burgers  equation.  Bewley  & 
Moin  (1994)  were  the  first  to  apply  the  suboptimal  control  procedure  to  a  turbulent 
flow  and  reported  about  17%  drag  reduction.  The  procedure  developed  by  Bewley 
&  Moin  (1994)  still  requires  velocity  information  inside  the  flow  in  order  to  solve 
the  adjoint  problem,  from  which  a  feedback  control  input  was  derived.  In  spite  of 
this  obvious  drawback,  however,  the  fact  that  a  control  theory  applied  to  a  turbulent 
flow  resulted  in  a  substantial  drag  reduction  is  encouraging,  since  their  control 
procedure  was  derived  rigorously  from  a  control  theory,  in  which  a  pre-determined 
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cost  functional  was  minimized.  Hill  (1993,  1994)  derived  a  control  input  as  a  function 
of  the  streamwise  wall  shear  only  by  modeling  the  near-wall  flow  with  a  spanwise 
velocity  growing  linearly,  and  normal  velocity  growing  quadratically,  with  normal 
distance  from  the  wall.  About  15%  drag  reduction  was  obtained  from  this  study. 
We,  however,  derive  a  simple  control  scheme  by  minimizing  cost  functionals  that 
are  related  to  the  streamwise  vortices,  which  have  been  found  to  be  responsible  for 
large  local  drag  in  turbulent  boundary  layers.  We  also  tried  to  minimize  drag  directly 
by  having  drag  itself  in  the  cost  functional,  but  it  was  not  successful  (see  §  3).  The 
objective  of  this  paper  is  to  demonstrate  that  a  wise  choice  of  the  cost  functional 
coupled  with  a  variation  of  the  formulation  can  lead  to  a  more  practical  control  law. 

We  present  how  to  choose  a  cost  functional  and  how  to  minimize  it  to  yield  simple 
feedback  control  laws  that  require  quantities  measurable  only  at  the  wall.  One  of 
the  laws  requires  spatial  information  on  the  wall  pressure  over  the  entire  wall  and 
the  other  requires  information,  also  over  the  entire  wall,  on  one  component  of  the 
wall  shear  stress.  We  then  derive  more  practical  control  schemes  that  only  require 
local  wall  pressure  or  local  surface  shear  stress  information,  and  show  that  they  work 
equally  well. 

2.  Suboptimal  procedure 

We  follow  a  similar  procedure  used  by  Choi  et  al.  (1993)  and  Bewley  &  Moin 
(1994).  The  problem  under  consideration  is  a  turbulent  channel  flow,  for  which  the 
governing  equations  are  the  Navier-Stokes  and  continuity  equations  with  the  no-shp 
boundary  condition; 

^  ^  (2.U 

dt  ^^dxj  dxi  Redxjdxj’ 

(2.2) 

OXj 

with 

Ui\y,  =  4>ix,z,t)5a,  (2-3) 

where  t  is  the  time,  xi,X2,X3  are  the  streamwise,  wall-normal,  and  spanwise  directions 
respectively,  Ui  are  the  corresponding  velocity  components,  p  is  the  pressure,  and 
Re  is  the  Reynolds  number,  and  the  control  input  is  the  wall-normal  velocity  at  the 
wall,  (f>. 

All  variables  are  non-dimensionalized  by  the  wall  shear  velocity,  m^,  and  the  channel 
half- width,  S.  We  also  use  interchangeably  x,y,z  for  xj  and  u,v,w  for  u,-.  Periodic 
boundary  conditions  are  imposed  in  the  streamwise  and  spanwise  directions.  The  flow 
rate  in  the  streamwise  direction  is  kept  constant,  and  the  drag  is  measured  by  the 
mean  pressure  gradient  necessary  to  maintain  the  constant  flow  rate. 

We  found  that  the  choice  of  the  cost  functional  to  be  minimized  is  critical  in 
the  performance  of  the  control.  Since  the  streamwise  vortices  have  been  known  to 
be  responsible  for  large  drag  in  turbulent  boundary  layers,  we  tried  to  choose  the 
cost  functional  that  is  directly  related  to  them.  This  is  based  on  our  conjecture  that 
a  suitable  manipulation  of  the  streamwise  vortices  would  lead  to  drag  reduction. 
We  carefully  selected  two  cost  functionals  based  on  our  observation  of  a  successful 
control  of  Choi  et  al.  (1994).  As  shown  in  figure  1,  Choi  et  al's  (1994)  blowing 
and  suction,  which  are  equal  and  opposite  to  the  wall-normal  velocity  component 
at  y'*'  =  10,  effectively  suppress  a  streamwise  vortex  by  counteracting  up-and-down 
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Roure  1.  Schematic  of  a  pressure  field  induced  by  a  control  based  on  y*"  =  10. 


Figure  2.  Correlation  between  a  pressure  field  of  a  no-control  case  and  a  pressure  field  modified 
by  a  control  based  on  information  at  y'*'  =  10.  A  line  indicating  no  change  of  the  pressure  field  is 
drawn  for  guidance. 


motion  induced  by  the  vortex.  This  blowing  and  suction  creates  locally  high  pressure 
in  the  near-wall  region  marked  with  and  low  pressure  in  the  region  marked  with 
’  in  figure  1.  A  crucial  aspect  of  the  present  analysis  is  the  observation  that  this 
blowing  and  suction  increases  the  pressure  gradient  in  the  spanwise  direction  under 
the  streamwise  vortex  near  the  wall.  To  demonstrate  this  behaviour,  we  examine 
computed  flow  fields.  Figure  2  shows  a  scatter  plot  between  two  pressure  gradient 
fields:  5p/3z|®  is  the  pressure  gradient  before  the  control  is  applied  and  dp/dz\f,  is 
the  pressure  gradient  at  the  same  location  after  the  control  of  Choi  et  al.  (1994) 
is  applied  for  one  time  step.  It  is  apparent  that  the  control  increased  the  pressure 
gradient  significantly.  For  the  uncontrolled  case,  based  on  the  wall-shear  velocity 
and  the  channel  half-width  is  about  110.  The  spectral  code  of  Kim,  Moin  &  Moser 
(1987)  is  used  for  all  the  computations  presented  here.  Simulations  are  carried  out 
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using  32  X  65  X  32  spectral  modes  (with  dealiasing  in  the  streamwise  and  spanwise 
directions)  with  the  computational  domain  of  4nd  x  25  x  4K5f3,  respectively. 

The  above  argument  suggests  that  we  should  seek  blowing  and  suction  that  increases 
the  pressure  gradient  in  the  spanwise  direction  near  the  wall  for  a  short  time  period 
(i.e,  in  the  suboptimal  sense)  in  order  to  achieve  a  similar  drag  reduction  to  that 
achieved  by  Choi  et  a/.’s  (1994)  control.  The  cost  functional  ^{(l>)  to  be  minimized  is 
then 

^  f  n  ....  1  /" 


dt  dS, 


where  the  integrations  are  over  the  wall  (S)  in  space  and  over  a  short  duration  in  time 
At,  which  typically  corresponds  to  the  time  step  used  in  the  numerical  computation, 
and  /  is  the  relative  price  of  the  control  since  the  first  term  on  the  right-hand  side 
represents  the  cost  of  the  actuation  (j).  Note  that  there  is  a  minus  sign  in  front  of 
the  second  term  since  we  want  to  maximize  the  pressure  gradient.  It  should  be  noted 
that  the  spanwise  pressure  gradient  at  the  wall  will  be  eventually  reduced  when  the 
strength  of  the  near-wall  streamwise  vortices  is  reduced  through  successful  control. 
Here,  blowing  and  suction  that  increase  the  spanwise  pressure  gradient  for  the  next 
step  are  sought  as  a  suboptimal  control. 

To  minimize  the  cost  functional,  we  first  define  the  differential  states  of  the  velocity 
and  pressure,  (0i,p),  using  a  Frechet  differential  (Finlayson  1972), 

(2.5) 


where 


am. 


/(0  +  e0) -/((/)) 


^  being  an  arbitrary  perturbation  field  to  4>. 

Next,  we  choose  the  Crank— Nicolson  scheme  for  the  linear  terms  and  a  Runge- 
Kutta  scheme  for  the  nonUnear  terms  to  yield  a  discretized  form  of  (2.1)  and  (2.2): 


At 

(2.8) 

2Re  dxjdXj 

2  dxi 

=  0, 

(2.9) 

dxj 

=  </>  5,1, 

(2.10) 

where  the  superscripts  n  + 1  and  n  denote  the  time  step,  and  R"  includes  the  nonlinear 
terms  and  the  explicit  parts  of  the  pressure  gradient  and  viscous  terms.  The  Frechet 
differential  of  (2.8)-(2.10)  yields  the  governing  equations  for  the  differential  states 


,  At 

'  2Re  dxjdxj  2  5x, 


(2.11) 


(2.12) 


with 
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er‘U  =  <?^/2.  (2.13) 

Note  that  $  =  0.  Hereinafter,  we  drop  the  superscript  n+1  and  all  variables 

are  understood  to  be  at  the  (n  +  l)th  time  step.  Note  that  there  is  no  contribution 
from  the  nonlinear  terms,  thus  making  the  equations  linear.  Generally,  the  suboptimal 
formulation  depends  on  the  time  advancement  scheme  used,  as  shown  here  (see  also 
Choi  et  al.  1993).  Dropping  the  nonlinear  terms  may  miss  important  flow  dynamics. 
However,  we  found  from  our  numerical  tests  with  the  full  nonlinear  terms  included 
that  the  contribution  from  the  nonlinear  terms  is  negligible  in  our  boundary  control 
with  short  optimization  interval  At;  it  turns  out  that  the  conservation  of  mass  due  to 
the  wall  actuation  dominates  the  near-wall  dynamics. 

Under  the  approximation  that  IRefAt  k^,  where  k  =  {kl  +  k^Y^^,  and  kx  and 
kz  denote  the  streamwise  and  spanwise  wavenumbers  in  the  x-  and  z-directions 
respectivelyf,  the  above  equations  have  the  following  solutions  in  the  semi-infinite 
domain  with  periodic  conditions  in  the  x-  and  z-directions  (see  the  Appendix): 


^iiy)  =  (exp[-(21?e/A0‘/"y]  -e-'=0 , 

(2.14) 

03(y)  =  y  1  {exp[-{2Re/Aty'^y]  -  e-*>) , 

(2.15) 

O  * 

11 

(2.16) 

(2.17) 

where  §j,p,  and  ^  are  the  Fourier  coefiicients  of  6j,p,  and  ^  respectively.  For  the 
channel  geometry,  we  originally  considered  both  walls  and  found  that  the  interaction 
between  two  walls  is  negligible  as  long  as  the  typical  wavelength  (~  2n/k)  associated 
with  near-wall  structures  is  much  smaller  than  the  channel  width. 

The  Frechet  differential  of  the  cost  functional  (2.4)  becomes 


^<l>  ^  AAt 


u: 


dt  dS 


-u[ 


dz 


The  Fourier  representation  of  the  above  equation  is 


=  4|-  _  «£ 


dp_ 

dz 


dz 


6tdS. 


(2.18) 


(2.19) 


where  the  hat  denotes  the  Fourier  coefficient,  and  the  superscript  •  denotes  the 
complex  conjugate.  From  (2.17),  dp/dz  can  be  expressed  in  terms  of 


dp 

dz 


2ikz 
:<p  • 


k  Ar 


Equation  (2.19)  then  reduces  to 


—4,  .14,4, 


Pw 


(2.20) 


(2.21) 


t  It  can  be  shown  that  ^  >  1.  Here  we  used  uAt/Ax  0(1)  and 

fcfmix  Re^^^y  where  ktj  is  the  wavenumber  corresponding  to  the  Kohnogorov  length  scale. 


250 


C.  Lee,  J.  Kim  and  H.  Choi 


dy  w 


Figure  3.  Correlation  between  a  spanwise  wall-shear  stress  of  a  no-control  case  and  a  spanwise 
wall-shear  stress  modified  by  a  control  based  on  y'^  =  10.  A  line  indicating  no  change  of  the 
spanwise  wall-shear  stress  field  is  drawn  for  guidance. 


For  an  arbitrary  $,  equation  (2.21)  should  be  satisfied,  yielding 


^<l> 


2kl 


(2.22) 


From  the  requirement  that  the  Frechet  differential  of  the  cost  functional  be  minimized, 
i.e.  ^ l&(f>  =  0,  the  optimum  ^  then  becomes 


i_r^h 

(p  —  C  ^  Pvv> 


(2.23) 


where  C  is  a  positive  scale  factor  that  determines  the  cost  of  the  actuation.  Equation 
(2.23)  indicates  that  the  optimum  wall  actuation  is  negatively  proportional  to  the  sec¬ 
ond  spanwise  derivative  of  the  wall  pressure,  with  the  high- wavenumber  components 
reduced  by  l/k. 

Another  wall  quantity  that  indicates  similar  changes  of  the  near-wall  dynamics  due 
to  the  altered  pressure  field  is  the  spanwise  shear  at  the  wall,  dw/dy.  Owing  to  the 
added  pressure  gradient  in  the  spanwise  direction  below  the  streamwise  vortex,  the 
spanwise  flow  near  the  wall  is  also  induced,  thus  increasing  the  spanwise  shear  stress 
at  the  wall  (see  figure  1).  The  scatter  plot  between  the  spanwise  shear  stress  at  the 
wall  from  two  different  fields,  one  of  which  is  from  an  unperturbed  channel  and  the 
other  with  the  control  based  on  y'^  =  10,  is  shown  in  figure  3.  Thus  another  choice 
for  the  cost  functional  to  be  minimized  is 


/{(!>)  = 


2AAt 


r  rf+At 

Jsl  ‘ 


4>^  dt  dS  - 


2AM 


JsJt  v^y/, 


dt  dS 


(2.24) 
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Rgure  4.  Correlation  between  a  pressure  field  of  a  no-control  case  and  a  pressure  field  modified 

by  a  control  based  on  equation  (2.23). 


Following  the  procedure  that  led  to  equation  (2.23)  yields  the  optimum  actuation: 


=  C 


i/cz  dw 
k  dy 


(2.25) 


Equation  (2.25)  indicates  that  the  optimum  wall  actuation  should  be  proportional  to 
the  spanwise  derivative  of  the  spanwise  shear  at  the  wall,  with  the  high-wavenumber 
components  reduced  by  1/fe.  Note  that  the  scale  factor  C  in  equations  (2.23)  and 
(2.25)  is  arbitrary.  As  mentioned  before,  the  suboptimal  formulation  depends  on  the 
time-advancement  scheme  used.  However,  if  we  used  a  different  implicit  scheme,  only 
the  resulting  constant  C  would  be  different.  This  does  not  cause  a  problem  since  C 
is  chosen  such  that  the  r.m.s.  value  of  the  wall  actuation  is  maintained  at  a  given 
value. 

In  the  following  simulations,  we  set  the  r.m.s.  value  of  to  be  equal  to  that  of  the 
wall-normal  velocity  at  y*"  =  10,  which  gives  the  same  r.m.s.  value  of  wall  actuations 
as  that  of  Choi  et  al.  (1994). 


3.  Results 

The  above  control  laws  (2.23)  and  (2.25)  are  tested  in  a  turbulent  channel  flow. 
The  pressure  gradient  at  the  wall  is  monitored  to  see  if  the  equation  (2.23)  control 
increases  the  pressure  gradient  when  the  control  law  is  applied.  It  behaves  as  expected, 
as  shown  in  figure  4.  The  spanwise  shear  stress  at  the  wall  also  increases  when  the 
equation  (2.25)  control  is  applied  (see  figure  5).  These  results  also  confirm  that  a 
suboptimal  procedure  without  the  nonlinear  terms  does  not  cause  an  error  for  our 
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figure  5.  Correlation  between  a  spanwise  wall-shear  stress  of  a  no-control  case  and  a  spanwise 
wall-shear  stress  modified  by  a  control  based  on  equation  (2.25). 


boundary  controls.  Since  these  control  laws  specifically  increase  the  wall  pressure 
gradient  or  spanwise  wall-shear  stress,  correlations  between  these  variables  before 
and  after  control  are  much  higher  than  those  corresponding  to  the  control  based  on 
y+  =  10  information.  This  suggests  that  the  controls  based  on  equations  (2.23)  and 
(2.25)  modify  the  flow  in  a  different  manner  from  the  control  based  on  information 
at  y"*"  =  10. 

Time  histories  of  the  mean  streamwise  wall-shear  stress  for  different  control  laws 
are  shown  in  figures  6  and  7.  As  control  begins,  an  immediate  drop  in  the  shear  stress 
is  observed  for  all  cases.  At  the  same  expense  (i.e.  the  same  r.ms.  value  of  (j)),  the 
control  based  on  (2.25),  which  reduces  drag  by  as  much  as  22%,  is  apparently  more 
effective  than  that  based  on  (2.23)  which  reduces  drag  by  16%.  This  indicates  that 
the  spanwise  wall-shear  stress  is  a  better  quantity  for  control  input. 

The  control  laws  presented  above,  however,  are  still  impractical  to  implement,  since 
they  are  expressed  in  terms  of  the  Fourier  coefficients  (i.e.  in  wavenumber  space), 
which  require  information  over  the  entire  spatial  domain.  Therefore,  the  inverse 
transforms  of  kj/k  and  ikz/k  are  sought  numerically  so  that  the  convolution  integral 
can  be  used  to  express  the  control  laws  in  physical  space.  The  discrete  representation 
of  each  control  law  then  becomesf: 


<l>{Xj,Zk)  —  C  ^  ^  ^  ,  Wyi^iPw{Xj+f,Zk+k') 
}'  k’ 


(3.1) 


t  Since  k^/k  and  ifej/fe  are  not  periodic  in  the  wavenumber  space,  a  high  resolution  was  used  to 
obtain  WJ^  and  IV)"  to  minimize  aliasing  error. 
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k{—  z/Az) 

0  1  2  3  4  5  6 

0  1.0000  -0.4427  0.0031  -0.0440  0.0044  -0.0138  0.0032 

7(=x/Ax)  1  0.0413  -0.0007  -0.0040  -0.0037  -0.0029  -0.0021  -0.0016 

2  -0.0110  0.0057  0.0019  0.0011  0.0002  0.0000  -0.0003 

Table  1.  The  weight  distribution  of  equation  (3.1),  The  weights  are  symmetric  between  j  and 

—j  and  k  and  —k.  The  weights  are  normalized  by  W^,  Bold  faced  weights  are  used  in  the  calculation 
of  figure  6.  Here,  =  40  and  Az"^  =  13  (Ax  =  3Az). 


k{=  z/Az) 

0  1  2  3  4  5  6 

0  0.0000  1.0000  -0.1039  0.2679  -0.0852  0.1419  -0.0671 

7(=x/Ax)  1  0.0086  0.0537  0.0503  0.0310  0.0340  0.0148  0.0237 

2  0.0001  -0.0104  0.0059  0.0051  0.0100  0.0074  0.0092 

Table  2.  The  weight  distribution  of  equation  (3.2),  The  weights  are  symmetric  between  j  and 
— 7  and  antisymmetric  between  k  and  —k.  The  weights  are  normalized  by  Bold  faced  weights 
are  used  in  the  calculation  of  figure  7.  Here,  Ax"^  =  40  and  Az"*"  =  13  (Ax  =  3Az). 


and 


f  U!  y 


{Xj+f,Zk+w), 

W 


(3.2) 


respectively.  The  subscripts,  j  and  k,  denote  the  discretizing  indices  in  the  x-  and 
z-directions  respectively.  The  weights,  Wf/^  and  WJ^,  are  given  in  tables  1  and  2.  Note 
that  the  weights  decay  rapidly  with  distance  from  the  point  of  interest,  suggesting  that 
the  optimum  actuation  can  be  obtained  by  a  local  weighted  average  of  the  pressure  or 
spanwise  shear  stress.  The  results  obtained  using  a  37-point  average  for  the  pressure 
and  an  11-point  average  for  the  spanwise  shear  stress  yield  about  the  same  drag 
reduction  as  that  obtained  from  full  integration  using  equations  (2.23)  and  (2.25)  (see 
figures  6  and  7).  The  weights  used  in  these  calculations  are  bold-faced  in  tables  1 
and  2.  It  is  remarkable  that  localized  information  can  produce  such  a  significant  drag 
reduction,  especially  for  the  control  with  the  spanwise  shear  stress. 

The  weight  distribution  WJ^  for  y  =  0  is  very  similar  to  the  one  found  in  the 
application  of  neural  networks  to  the  same  turbulent  flow  (Lee  et  al.  1997),  in  which 
only  a  single  strip  of  the  spanwise  shear  information  was  used  in  the  network.  The 
control  law  of  Lee  et  al., 


-  _  ikz  3w 

\kz\  dy 


(3.3) 


produces  almost  the  same  distribution  of  blowing  and  suction  as  the  present 
results.  Note  that  blowing  and  suction  from  equation  (2.25)  are  nearly  the  same 
as  those  from  equation  (3.3)  because  the  near-wall  structures  have  relatively  slow 
variation  in  the  streamwise  direction  (the  kx  =  0  component  is  dominant).  This  blow¬ 
ing  and  suction  distribution  is  also  very  similar  to  the  one  based  on  y*"  =  10  data 
(see  Lee  et  al.  1997). 

In  figure  8,  contours  of  the  streamwise  vorticity  in  a  cross-flow  plane  for  the  control 
with  equation  (2.25)  are  compared  with  a  no-control  case.  Significant  reduction  in 
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Figure  6  Mean  streamwise  wall-shear  stress  for  various  control  laws  based  on  the  wall  pressure 
compped  to  the  no-control  case:  thick  solid  line,  the  control  law  expressed  in  Fourier  space,’ 
equation  (2.23);  thin  solid  line,  the  control  law  expressed  in  physical  space,  equation  (3  1)  with  the 
integration  radius  of  6Az.  Ar**"  0,2. 


Rgum  7.  Mean  streamwise  wall-shear  stress  for  various  control  laws  based  on  the  spanwise 
wall-shear  stress  compared  to  the  no-control  case:  thick  solid  line,  the  control  law  expressed  in 
®4uation  (2.25);  thin  solid  line,  the  control  law  expressed  in  physical  space,  equation 
(3.2)  with  1 1  points  in  the  spanwise  direction  only.  At'*'  ~  0.2. 


the  strength  of  the  streamwise  vortices  is  evident.  This  again  confirms  the  notion 
that  a  successful  manipulation  of  the  near-wall  streamwise  vortices  can  lead  to  drag 
reduction.  The  mean  streamwise  velocity  and  root-mean-square  velocity  near  the  wall 
are  shown  in  figures  9  and  10,  respectively,  and  compared  with  the  no-control  case. 
The  trends  are  very  similar  to  Choi  et  al.  (1994)  and  Lee  et  al.  (1997). 

Finally,  we  mention  that  we  tried  to  find  a  feedback  control  scheme  minimizing  the 
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Rgure  8.  Contours  of  the  streamwisc  vorticity  in  a  cross-flow  plane;  (a)  no  control;  (b)  control 
based  on  wall-shear  stress  (equation  (2.25)).  The  contour  level  increment  is  the  same  for  both  figures. 
Negative  contours  arc  dashed. 


figure  9.  The  mean  stream  wise  velocity  normalized  by  u^:  Thick  solid  line,  control  law  (2.25); 
dashed  line,  no  control.  For  the  controlled  case,  the  velocity  is  normalized  by  the  controlled  Uj. 
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Figure  10.  The  root-mean-square  velocity  normalized  by  of  the  uncontrolled  flow:  thick  solid 
line,  control  law  (2.25);  dashed  line,  no  control. 


following  cost  functional: 


/(0)  = 


2AAt 


dt  dS  ”|- 


1  f 


dt  dS, 


(3.4) 


with  m  =  1  or  2.  This  cost  functional  is  the  most  natural  choice  since  it  contains 
the  quantity  directly  related  to  drag.  We  followed  the  same  procedure  to  derive  the 
following  control  schemes: 

.^  =  0,  (3.5) 


for  m  =  1,  and 


ikx  du 
k  dy 


W 


(3.6) 


for  m  =  2.  Equation  (3.5)  obviously  does  not  reduce  drag  and  equation  (3.6)  increased 
drag  in  our  numerical  simulations.  This  failure  appears  to  be  due  to  the  neglect  of 
the  nonlinear  terms  in  our  formulation,  since  the  numerical  solution  of  the  optimum 
wall  actuations  with  the  full  nonlinear  terms  gave  different  results  (Bewley  &  Moin 
1994).  This  suggests  that  even  in  a  short  time  interval  the  nonlinear  terms  should 
be  included  in  the  suboptimal  formulation  when  drag  itself  is  chosen  as  the  cost 
functional.  Manipulation  of  the  streamwise  vortices  by  wall  actuation,  however,  can 
be  accompUshed  through  a  linear  process  as  shown  in  the  previous  section.  It  appears 
that  having  a  term  that  is  directly  related  to  near-wall  streamwise  vortices  in  the  cost 
functional  indirectly  includes  the  nbnUnear  effect.  This  is  probably  related  to  the  fact 
that  near-wall  streamwise  vortices  are  self-sustained  by  a  nonUnear  process  (Hamilton, 
Kim  &  Waleffe  1995).  If  the  nonlinear  terms  are  included,  of  course,  it  is  impossible 
to  derive  the  feedback  control  law  in  closed  form  without  approximation.  Hill  (1993), 
on  the  other  hand,  considered  the  nonhnear  term  by  modelUng  the  near-wall  flow 
using  the  Taylor  series  expansion  and  presented  a  feedback  control  law  in  closed 
form,  which  resulted  in  about  15%  drag  reduction. 
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4.  Summary 

We  have  obtained  two  simple  feedback  control  laws  for  drag  reduction  by  applying 
a  suboptimal  control  procedure  to  a  turbulent  flow.  This  was  possible  since  we  selected 
two  cost  functionals  guided  by  the  successful  control  based  on  quantities  monitored 
at  y*"  =  10.  The  present  control  laws  perform  very  well,  resulting  in  substantial 
drag  reductions  when  applied  to  a  turbulent  channel  flow.  More  convenient  control 
schemes  requiring  only  local  information  were  also  derived,  and  were  shown  to  work 
equally  well.  They  require  quantities  measurable  only  at  the  surface  and  thus  should 
be  easier  to  implement  in  practice.  The  present  results  further  substantiate  the  notion 
that  a  successful  manipulation  of  the  near-wall  streamwise  vortices  is  the  key  to 
boundary-layer  control  for  drag  reduction. 

We  are  grateful  to  Dr  Gary  Coleman  and  Mr  Hong  Zhang  for  useful  discussions. 
This  work  was  supported  by  AFOSR  grant  F49620-95-10263  (Program  Managers, 
Drs  James  M.  McMichael  and  Marc  N.  Glauser).  Computer  time  was  provided  by 
the  NAS  Program  of  NASA  Ames  Research  Center,  the  San  Diego  Supercomputer 
Center,  and  the  DoD  Eflgh  Performance  Computing  Center.  H.C.  acknowledges  the 
support  by  KOSEF  through  Turbo  and  Power  Machinery  Research  Center. 


Appendix.  Derivation  of  the  solution  (2.14)-(2.17) 

Equations  for  the  differential  states  (ObP),  (2.11),  (2.12)  can 
of  the  Fourier  coefiicients. 

be  rewritten  in  terms 

(Al) 

(A  2) 

(A3) 

ikx0{  +  H"  =  Oj 

dj; 

(A4) 

with 

k0)  =  $5a,  koo)  =  0,  (A  5) 

where  and  p  are  defined  as  follows : 

=  (A  6) 

k,  k, 

=  (A  7) 

k,  k. 

The  operation  ifc*-  (A  1)  -l-d/dy(A2)  -fikz-  (A3)  together  with  equation  (A 4)  yields 


0. 


which  has  a  non-growing  solution. 


(A8) 

(A  9) 
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with  a  wall  value,  p„,  which  will  be  determined  later.  Using  (A  9),  we  can  find  solutions 
to  equations  (A  1),  (A 2),  (A3), 

6iiy)  =  y  (exp[-(fc^  +  2Re/At)^^^y]  -  e“‘^) ,  (A  10) 

^sCy)  =  y (exp[-(/c^  +  2Re/Aty'^y]  - e”*'^) ,  (All) 

^2(y)  =($-  jkp„^  exp[-(fe'  +  2Re/Aty/^yJ  +  yfc^.e-*=^.  (A  12) 
With  these,  equation  (A  4)  reduces  to 

(“  (*'  +  2RelAtyi-^y\  =  0.  (A  13) 

Since  2Rel At  »  fc^,  equation  (A  13)  yields 

^  =  ^1.  (A  14) 

With  this,  equations  (A  10),  (All),  (A  12),  (A9)  become  equations  (2.14),  (2.15),  (2.16), 
(2.17),  respectively. 
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A  systems  theory  framework  is  presented  for  the  linear  stabilization  of  two-dimen¬ 
sional  laminar  plane  Poiseuille  flow.  The  governing  linearized  Navier-Stokes  equations 
are  converted  to  control-theoretic  models  using  a  numerical  discretization  scheme. 
Fluid  system  poles,  which  are  closely  related  to  Orr-Sommerfeld  eigenvalues,  and 
fluid  system  zeros  are  computed  using  the  control-theoretic  models.  It  is  shown 
that  the  location  of  system  zeros,  in  addition  to  the  well-studied  system  eigenvalues, 
are  important  in  linear  stability  control.  The  location  of  system  zeros  determines  the 
effect  of  feedback  control  on  both  stable  and  unstable  eigenvalues.  In  addition,  system 
zeros  can  be  used  to  determine  sensor  locations  that  lead  to  simple  feedback  control 
schemes.  Feedback  controllers  are  designed  that  make  a  new  fluid-actuator-sensor— 
controller  system  linearly  stable.  Feedback  control  is  shown  to  be  robust  to  a  wide 
range  of  Reynolds  numbers.  The  systems  theory  concepts  of  modal  controllability  and 
observability  are  used  to  show  that  feedback  control  can  lead  to  short  periods  of  high- 
amplitude  transients  that  are  unseen  at  the  output.  These  transients  may  invalidate  the 
linear  model,  stimulate  nonlinear  effects,  and/or  form  a  path  of ‘bypass’  transition  in  a 
controlled  system.  Numerical  simulations  are  presented  to  validate  the  stabilization  of 
both  single-wavenumber  and  multiple-wavenumber  instabilities.  Finally,  it  is  shown 
that  a  controller  designed  upon  linear  theory  also  has  a  strong  stabilizing  effect  on 
two-dimensional  finite-amplitude  disturbances.  As  a  result,  secondary  instabilities  due 
to  infinitesimal  three-dimensional  disturbances  in  the  presence  of  a  finite-amplitude 
two-dimensional  disturbance  cease  to  exist. 


1.  Introduction 

A  basic  problem  in  fluid  dynamics  is  the  theoretical  understanding  of  how  insta- 
bihty  in  laminar  shear  flow  leads  to  transition  to  turbulence.  Since  laminar  flow 
is  preferred  in  many  applications,  the  suppression  of  fluid  instabilities  to  maintain 
laminar  flow  would  be  very  useful.  Towards  this  goal,  active  boundary  layer  control 
of  instability  has  been  proposed.  The  nonlinear  aspects  of  the  transition  process  are 
still  not  clearly  known.  It  has  been  shown  (Orszag  &  Patera  1983)  that  some  shear 
flows  may  sustain  two-dimensional  finite-amplitude  instabilities  that  cause  infinites¬ 
imal  three-dimensional  disturbances  to  be  highly  unstable.  This  two-dimensional 
primary  instability /three-dimensional  secondary  instability  process  may  be  one  non¬ 
linear  mechanism  that  leads  to  transition.  Conversely,  the  behaviour  of  infinitesimal 

t  Present  address:  Jet  Propulsion  Laboratory,  4800  Oak  Grove  Drive,  Pasadena,  CA  91109,  USA. 
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(linear)  disturbances  in  laminar  shear  flow  is  well  understood  (Orr  1907;  Sommerfeld 
1908;  Drazin  &  Reid  1981).  Therefore,  linear  analysis  provides  a  natural  starting 
point  to  begin  to  develop  active  control  schemes  that  may  eventually  lead  to  full 
transition  control.  It  is  the  purpose  of  this  paper  to  develop  feedback  controllers 
based  on  linear  theory  that  stabilize  two-dimensional  plane  Poiseuille  flow  to  in¬ 
finitesimal  disturbances.  In  addition,  it  will  be  shown  that  a  controller  designed 
upon  linear  theory  has  a  strong  stabilizing  effect  on  two-dimensional  finite-amplitude 
instabilities.  As  a  result,  secondary  three-dimensional  instabilities  as  described  in 
Orszag  &  Patera  (1983)  cease  to  exist  in  such  a  system. 

Most  prior  work  in  the  area  of  laminar  flow  linear  instability  suppression  has 
concentrated  on  the  wave  superposition  approach.  A  nice  survey  of  past  work  is 
given  in  Joslin,  Erlebacher  &  Hussaini  (1994).  The  basic  idea  is  that  boundary  layer 
instabilities  appear  as  a  combination  of  many  sinusoidally  growing  waves  of  certain 
frequencies,  phases,  and  amplitudes.  If  these  wave  parameters  are  known,  or  if  they 
can  be  determined,  a  second  wave  may  be  stimulated  in  the  flow  that  is  exactly  out 
of  phase  with  the  instability  wave.  In  this  way,  the  two  waves  may  destructively 
interfere  and  flow  may  be  stabilized.  Disturbance  of  the  flow  field  to  cause  a 
wave  motion  to  appear  has  been  experimentally  demonstrated  using  several  methods 
such  as  vibrating  ribbons  (Schubauer  &  Skramstad  1947),  vibrating  wires  (Milling 
1981),  and  heating  elements  (Nosenchuck  1982).  In  addition,  several  authors  have 
numerically  obtained  wave  superposition  results  based  on  Navier— Stokes  simulations 
(Beringen  1984;  Metcalfe  et  al.  1986).  Most  of  the  results,  however,  report  incomplete 
destruction  of  instability  waves.  Joslin  et  al.  (1994)  explain  that  wave  cancellation 
is  very  sensitive  to  the  wave  parameters  and  postulate  that  incomplete  destruction 
reported  in  past  studies  was  due  to  improper  phase  or  amplitude  properties  of  the 
cancelling  wave. 

In  addition.  Bower,  Kegelman  &  Pal  (1987)  considered  the  Orr-Sommerfeld  equa¬ 
tion  in  designing  an  input  to  channel  flow  that  may  counteract  the  effects  of  a 
disturbance  that  excites  instabilities.  They  showed  that  if  an  oscillating  flat  pulse 
function  at  the  lower  boundary  of  a  channel  excited  inherent  instabiUties,  a  second 
oscillating  flat  pulse  function  could  be  constructed  downstream  that  negated  excite¬ 
ment  of  the  instabilities.  Much  like  the  wave  superposition  approach,  their  aim  was 
not  to  stabilize  the  underlying  dynamics  of  the  problem,  but  rather  to  ‘cancel’  the 
effects  of  a  specific  disturbance. 

Unlike  past  work,  our  aim  is  to  use  systems  theory  to  construct  a  combination  fluid- 
actuator-sensor-controller  system  that  is  inherently  stable.  In  essence,  the  approach 
changes  the  philosophy  of  the  problem  from  thinking  about  how  inputs  can  mitigate 
an  inherently  unstable  system  to  thinking  about  how  sensors  and  actuators  can  be 
added  to  form  an  entirely  new  stable  system.  Recently,  mathematical  control  systems 
theory  has  begun  to  be  applied  to  fluid  systems  (Burns  &  Ou  1994;  Gunzburger, 
Hou  &  Suobodny  1992;  Choi,  Moin  &  Kim  1993).  A  control  theory  approach  for 
laminar  flow  linear  instability  suppression  will  be  shown  to  have  several  advantages 
over  the  traditional  wave  cancellation  approach.  It  will  eliminate  the  need  to  explicitly 
measure  phase  and  frequency  of  instabilities.  Also,  it  will  provide  a  framework  to 
select  sensor  locations  in  order  to  have  the  least  complex  controller.  Further,  feedback 
control  will  be  shown  to  be  extremely  robust  to  changing  Reynolds  numbers  given 
proper  sensor  location.  In  addition,  linear  feedback  controllers  will  be  shown  to  have 
a  strong  stabilizing  effect  on  two-dimensional  finite-amplitude  disturbances. 

This  paper  is  organized  as  follows.  In  §2,  we  formulate  the  linear  channel  flow 
problem  using  the  linearized  Navier-Stokes  equations.  Boundary  input  in  the  form 
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of  blowing/suction  and  boundary  output  in  terms  of  shear  are  represented  within  the 
equations.  In  §3,  we  introduce  the  concepts  of  zeros  and  poles  of  a  system,  as  well  as 
control-theoretic  models.  The  governing  partial-differential  equation  for  the  system 
is  converted  into  a  set  of  first-order,  ordinary  differential  equations  via  a  Galerkin 
method.  These  first-order,  ordinary  differential  equations  are  then  converted  into  a 
control-theoretic  model.  Section  4  describes  the  infinite-dimensional  nature  of  the 
channel  flow  system  and  how  it  affects  the  selection  of  actuation.  Section  5  describes 
the  numerical  model  and  the  verification  of  the  calculated  poles  and  zeros.  Section  6 
introduces  feedback  control  design.  It  is  shown  that  judicious  sensor  placement,  based 
on  zero  locations,  can  lead  to  simple  control  schemes.  Furthermore,  the  control  system 
is  extremely  robust  to  change  in  Reynolds  number.  Section  7  explores  unobservable 
modal  reinforcement  as  a  possible  path  of  ‘bypass’  transition  in  a  controlled  system. 
It  then  shows  how  a  particular  control  model,  called  the  modal-canonical  state 
space  model,  may  be  used  to  assess  observability  of  modal  reinforcement.  Section  8 
demonstrates  multiple-wavenumber  instability  control.  Section  9  demonstrates  that  a 
linear  controller  has  a  strong  stabilizing  effect  on  two-dimensional,  finite-amplitude 
instabilities.  As  a  result,  secondary  three-dimensional  instabilities  as  described  by 
Orszag  &  Patera  (1983)  cease  to  exist.  Finally,  §10  outlines  conclusions. 


2.  Problem  formulation 

2.1.  Dynamic  equations 

We  consider  two-dimensional  plane  Poiseuille  flow  between  two  parallel  stationary 
plates.  Let  the  channel  be  of  finite  length  and  finite  height,  with  the  centreline  at 
zero.  The  flow  in  the  channel  is  described  by  the  unsteady  nonlinear  incompressible 
Navier-Stokes  equations.  In  order  to  study  the  linear  stabiUty  of  the  system,  we 
follow  the  standard  procedure.  Consider  small  perturbations  in  the  velocities  of 
u{x,y,t)  in  the  horizontal  direction,  v{x,y,t)  in  the  vertical  direction,  and  p(x,y,t)  in 
the  pressure  field.  Let  the  primary  flow  be  represented  by  U(y)  with  Uc  being  the 
centreline  velocity.  The  linearized  incompressible  Navier-Stokes  equations  may  be 
formed  by  substituting  the  primary  flow  and  small  perturbations  into  the  nonlinear 
incompressible  Navier-Stokes  equations  and  disregarding  the  second-order  terms 
involving  the  perturbations. 


ot  OX  oy  Re 

du(x,y,t)  dv{x,y,t)  _• 
dx  dy 


(2.1) 

(2.2) 

(2.3) 


where  the  flow  variables  are  non-dimensionalized  by  the  channel  half-height,  H,  and 
centreline  velocity,  Uc.  Re  is  the  Reynolds  number  defined  as  (UcH/v)  where  v  is  the 
kinematic  viscosity.  By  introducing  a  ‘stream  function’,  y>{x,  y,  t), 


dy 


(2.4) 


vix,y,t)  = 


dy}{x,y,t) 

dx 


(2.5) 


and 
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(2.1)-(2.3)  may  be  combined  into  a  single  equation, 


dtdx^  dtdy^  dx^ 


(2.6) 


Assume  periodic  boundary  conditions  in  the  streamwise  (x)  direction.  For  channel 
flow,  with  rigid  plates  at  y  =  —  1  and  y  =  1,  the  no-slip  boundary  conditions  become 


ip(x,y  = -l,t)  =  0,  (2.7) 

^(x,y  =  -l,t)  =  0,  (2.8) 

dy 

tpix,y=Ut)  =  0,  (2.9) 

^(x,y  =  l,t)  =  0.  (2.10) 

dy 

With  an  initial  condition 

t/)(x,  y,  r  =  0)  =  g(x,  y)  (2. 1 1) 


the  boundary  value  problem  is  completely  formed.  Existence  and  uniqueness  of  solu¬ 
tions  for  the  linearized  Navier-Stokes  equations  have  been  studied  in  Ladyzhenskaya. 
(1969),  Kreiss  &  Lorenz  (1989),  and  Temam  (1984).  Equations  (2.6)-(2.11)  represent 
the  starting  point  for  construction  of  a  feedback  control  system.  These  equations 
neither  include  any  control  terms  nor  do  they  describe  any  sensing  of  flow  field 
variables. 


2.2.  Boundary  input 

We  consider  the  case  of  blowing/suction  at  the  lower  wall  of  the  channel.  The  bound¬ 
ary  conditions  are  now  modified  from  before  to  include  boundary  input,  represented 


as  the  known  separable  function  q{t)w{x)f{y), 

w(x,y  =  -1.  t)  =  q(t)w(.x)f{y  =  -1),  (2.12) 

|^(x,y  =  -l.t)  =  q{t)w{x)^-~ — ^  =  0,  (2.13) 

i/)(x,y  =  l,r)  =  0,  (214) 

(2.15) 

Note  that  these  conditions  constrain  the  function  /(y)  such  that 

/(y  =  -1)  ^  0,  (2.16) 

=  0.  (2-17) 

dy 

f{y  =  1)  =  0,  (2.18) 

=  0.  (2.19) 

Many  functions  may  be  equally  appropriate.  One  such  function  is 

f(y)  =  \y*  +  \y^-y^-ly  +  ^-  (2-20) 
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In  order  to  relate  boundary  conditions  on  ip  to  blowing/suction  in  the  wall-normal 
direction,  we  use  (2.5)  to  relate  v(x,y,t)  and  ip(x,y,t).  Then  (2.12)  becomes 

«(^.  y  =  -1. 0  =  (2.21) 

Note  that  v{x,  y,  t)  is  related  to  the  derivative  of  w(x).  The  homogeneous  equation 
(2.6)  and  the  inhomogeneous  boundary  condition  (2.12)  can  be  converted  into  an 
inhomogeneous  equation  with  homogeneous  boundary  conditions  by  introducing 

<^(x,  y,  t)  =  ip(x,  y,  t)  -  q(t)fiy)w(x).  (2.22) 


Then  by  substituting  into  (2.6),  we  obtain 
5  5V.  d  5V  ru  ru,.,  d  5V  ,  d^Uiy)d<l>  ,  1  5V 

dt  5x2  +  St  5y2  ^^^5x2  ^^^5x  5y2  dy2  dx  Re  dx*  Redx^  5y2 


5,(0  _  ,(0^C/(y)/(y) 


1  5V 

Re  dy* 


dt  5x2 


dt  ^  '  dy 
5w(x)  d^l7(y) 


1 


5x2 
5‘*w(x) 


1 


dx  5y2 

5^w(x)  d^f{y) 


dx 


1 


dy2 


/W  +  r;9(') 


f(y) 


Sy2  +rJ9(«)»’W  gy. 


The  boundary  conditions  in  terms  of  (j)  are  now 

<^(y  =  -l)  =  0, 
5</.(y  =  -!)_, 
dy 

<l>(y  =1)  =  0, 
5<^(y=l) 


dy 


=  0. 


(2.23) 

(2.24) 

(2.25) 

(2.26) 
(2.27) 


2.3.  Boundary  output 

We  use  the  streamwise  component  of  shear  at  a  single  boundary  point,  z(xi,  y  =  —  1,  t), 
as  our  boundary  output,  which  is  given  by 

du{xi,y  =  -1,0 


z(Xi,y  =  -l,0 


dy 


By  expressing  u(xi,y  =  — 1,0  in  terms  of  the  stream  function  (2.4), 

52ip(xj,y  =  -1,0 


z(xi,y  =  -l,t)  = 


and  by  observing  (2.22) 

d^ip(xi,y  =  -l,t) 


z(Xi,y  =  -1,0  = 


5y2 


5y2 

5V(x,-,y  =  -1,0 

5^2 


(2.28) 


(2.29) 


5y2 


(2.30) 


3.  Zeros,  eigenvalues,  and  control-theoretic  models 

Linear  stability  analysis  of  (2.6)  (Drazin  &  Reid  1981)  shows  that  the  system  is 
linearly  unstable  for  a  range  of  Reynolds  numbers.  The  goal  of  this  paper  is  to 
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stabilize  the  system  using  control  theory.  To  do  this,  we  first  transform  the  governing 
equations  into  special  control-theoretic  models. 


3.1.  System  transfer  function 

We  have  defined  a  single-input/single-output  (SISO)  system  in  the  sense  that  only 
one  scalar  function,  q{t),  defines  the  entire  input  and  one  scalar  function,  z{t),  defines 
the  output.  A  common  form  of  control  model  for  a  finite-dimensional  SISO  system 
is  the  transfer  function  model.  The  transfer  function,  H{s),  is  defined  as  the  Laplace 
transform  of  the  output,  z(t),  divided  by  the  Laplace  transform  of  the  input,  q{t), 
where  zero  initial  conditions  are  assumed.  Then 


A  nz(t)]  Zis) 
nqit)]  QisY 


(3.1) 


For  finite-dimensional  systems,  Z(s)  and  (2(s)  take  the  form  of  polynomials  in  the 
complex  variable  s.  These  polynomials  may  be  factored  to  yield  an  equivalent 
representation. 


His) 


A  Z(s) 

Qis) 


J 


Uis-Cj) 

7=1 


/ 


n(s-Pi) 


(3.2) 


In  this  form,  pi . .  .p/  are  the  poles  of  the  system.  The  poles  of  any  system  are  dependent 
solely  on  the  physics  of  the  underlying  system,  independent  of  any  particular  input  or 
output.  Unstable  modes  of  the  system  appear  as  poles  whose  real  part  is  greater  than 
zero.  As  will  be  seen  in  later  sections,  the  poles  are  closely  related  to  the  eigenvalues 
of  the  Orr-Sommerfeld  equation.  The  values  Ci ...  0  are  the  zeros  of  the  system.  They 
are  heavily  dependent  on  which  particular  inputs  and  outputs  are  used  on  the  system. 
As  will  be  seen  later,  the  position  of  these  zeros  will  dictate  sensor  locations  and  will 
reveal  the  effect  of  feedback  control  on  the  eigenvalues.  A  graphical  representation  of 
the  transfer  function  can  be  produced  by  plotting  the  poles  and  zeros  in  the  complex 
s-space. 


3.2.  State-variable  model 

Much  of  modem  control  theory  is  based  on  the  state-variable  representation  of  a 
dynamic  system.  This  representation  relies  on  the  basic  fact  that  the  motion  of  any 
fiiiite-dimensional  dynamic  system  may  be  expressed  as  a  set  of  first-order  ordinary 
differential  equations.  As  a  simple  example  of  a  state  variable  model  (Franklin,  Powell 
&  Emami-Naeini  1988),  Newton’s  law  for  a  constant  single  mass,  M,  moving  in  one 
dimension,  x,  under  a  force,  F{t),  is 

M^  =  F{t).  (3.3) 


If  we  define  one  state  variable  as  the  position  xi  =  x(t)  and  the  other  state  variable 
as  the  velocity  X2  =  dx(0/dt,  (3.3)  can  be  written  as 


dx2  _  F{t) 
dt  M  ■ 


(3.4) 

(3.5) 
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Furthermore,  the  first-order  linear  ordinary  differential  equations  can  be  expressed 
using  matrix  notation 


dxi 

IT 

’or 

■  xi  ■ 

4_ 

■  0  ■ 

dX2 

0  0 

X2 

1 

1 

.  "dT  . 

If  we  take  output  as  position,  xi. 


(3.6) 


(3.7) 


z  =  [  1  0  ] 


Xi 

X2 


(3.8) 


or 

z  —  C '  X  (3.9) 

The  matrix  A,  and  the  vectors  B  and  C  are  called  the  state  space  matrices  of  the 
single-input,  single-output  (SISO)  system.  More  specifically,  the  A  matrix  is  referred 
to  as  the  dynamic  matrix  of  the  system.  It  can  be  shown  that  the  poles  of  the  system 
are  simply  the  eigenvalues  of  the  A  matrix.  In  the  more  general  case  of  multiple-input, 
multiple-output  (MIMO)  systems,  B  and  C  are  matrices.  For  generality,  we  may  add 
another  term,  Dq,  to  the  output  equation  to  account  for  systems  in  which  there  is 
direct  feedthrough  from  the  input  to  the  output.  Then, 

z  =  C-x  +  Dq  (3.10) 


In  the  case  of  no  direct  feedthrough,  D  —  0,  the  state  space  and  transfer  function 
models  are  related  as 


(3.11) 


3.3.  State-space  formulation  for  channel  system 

We  will  convert  our  problem  into  a  set  of  first-order  ordinary  differential  equations 
and  then  form  a  state-space  model  from  these  equations.  The  state-space  model  can 
then  be  represented  with  transfer  function  poles  and  zeros.  We  will  proceed  in  this 
way  for  two  reasons.  First,  our  system  lends  itself  to  decomposition  into  first-order 
ordinary-differential  equations  by  use  of  a  Galerkin  method.  More  importantly, 
however,  the  state-space  model  lends  itself  to  extremely  elegant  ways  to  control 
eigenvalues  of  a  system  in  well  prescribed  ways.  It  should  be  noted  that  unlike  the 
single-mass  example  given  above,  the  channel  system  requires  an  infinite  number  of 
ordinary  differential  equations  to  describe  its  motion.  This  is  known  as  an  infinite¬ 
dimensional  system.  As  a  result,  any  finite  number  of  ordinary-differential  equations 
used  in  a  state-space  model  will  not  completely  describe  the  system.  The  difficulties 
associated  with  such  a  system  are  taken  up  after  a  discussion  of  the  Galerkin  method 
used  to  obtain  the  ordinary-differential  equations. 


3.3.1.  First-order  system 

We  use  a  standard  Galerkin  procedure  to  convert  the  governing  partial  differ¬ 
ential  equation  (2.23)  into  a  system  of  ordinary-differential  equations.  Assume  an 


(3.12) 
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approximate  solution  of  (2.23)  as 

N  M 

n=-yv  m=0 

By  using  Fourier  functions,  for  P„{x)  and  basis  functions  constructed  from 

Chebyshev  polynomials  for  r„(y)  (Joshi  1996),  we  obtain  a  first-order  system.  Define 
inner  products  in  the  stream  wise  (x)  and  normal  (y)  directions  respectively: 

“j  rr«/2  ^  — 

=  f  /  e‘'^*e-™“<'Mx  =  oo  =  — ,  (3.13) 

L  J-i/2  L, 

and 

[PM  r„{y)]y  ^  jT*  (3.14) 

where  L  is  the  non-dimensional  length  of  the  finite-length  channel.  Applying  the 
orthogonality  of  the  basis  functions  in  x,  we  obtain  a  system  of  first-order  ordinary 
differential  equations; 

M  p  M  A  y  V  Af  M 

-  E  A + E  -  E  -  E 

m==0  m=0  m=0  m=0 

M  1  ^  1  ^ 

+  E«/m(t)««ojS^  +  ^  -  2^  ^aUt)l^c^Plk 

m=0  m=0  m=0 

+  ^  E  l  =  -N...N,k  =  0...M,  (3.15) 

m=0 
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(3.21) 


Jy 


3.3.2.  Transformation  to  state-space  form 
We  may  visualize  (3.15)  in  block  matrix  form  as 

^l=-N 

r  M.n  0  •  •  •  0  1  ' 

0  0  0 

0  0  0 


dfl/=-jv+i 

dt 


0 

0 

0 

Mn 

J 

1  dai^N 

L  dt 

K-n 

0 

0  1 

ai=-N 

0 

K-n+i 

0 

0 

a/=_jv+i 

0 

0 

0 

* 

0 

0 

0 

Kn  . 

+  [  r,  Y2' 


Q{t) 

dqjt) 

L  dt 


(3.22) 


where  a;=p  denotes  a  column  vector  of  all  aik  coefficients  whose  first  index,  I,  is  the 
value  p.  In  compact  notation, 


=  Y2] 


qit) 

dq{t) 


L  dt  J 

Assuming  non-singularity  of  the  M  matrix,  we  may  invert  to  obtain 


^  =  M-^Ka-i-  [  M-^Yi  M-^Y2  ] 


qit) 
dg(0 
dt  J 


(3.23) 


(3.24) 


where  a,  K,  Ti  and  Y2  are  all  complex.  By  expanding  in  terms  of  real  and  imaginary 
parts,  we  obtain 

^  4-  i^  =  M-^KrOr  M-^KiOr  -1-  iM-^KRtti  -  M-^Kja, 
at  at 

-\-M-%Rq{t)  4-  iM-^Yuqit)  4-  4-  iM-'rz/^.  (3.25) 

where  the  subscript  Rot  I  indicates  real  or  imaginary  parts.  Define 


-  A 

P  = 


«/ 

qit) 


(3.26) 


Then  the  state  space  system  is 


dt 

We  write  (3.27)  as 


■  m-^Kr 

s 

7 

S 

1 

■  M-^Y2R  ' 

= 

M-'Ki 

M-^Kr 

P  + 

M-^Y2I 

0 

0 

0 

1 

dqjt) 

dt 


^.Ap  +  B^&Ap  +  Bm 


(3.27) 


(3.28) 
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where  the  input  to  the  system  is  d.q{t)/At.  Note  since  q{t)  in  (2.12)  is  related  to 
physical  blowing/suction,  it  is  important  to  keep  in  mind  that  the  input  to  this  model 
is  the  derivative  of  q{t).  In  physical  terms,  any  input  derived  from  this  model  will 
be  in  terms  of  dq{t)/dt  and  therefore  must  be  integrated  before  it  can  be  applied  as 
physical  blowing/suction. 


3.3.3.  Measurement  equations 

By  observing  (2.30)  and  using  the  assumed  solution  form  (3.12) 

/  1  .N  d^xpixuy  = -l,t)  = 

z(x„  y - 1,  t)  - - - ^  +  q{t) - - wW 


n=—N  m=0 


(3.29) 


where  p(t,  x,-,  y  =  —  1),  the  residual  component  of  streamwise  shear  due  to  the  truncated 
expansion,  is  assumed  to  be  negligible  compared  to  the  two  other  terms.  The  second 
term  on  the  right-hand  side  of  (3.29)  is  made  up  of  the  known  input  terms  (2.12). 
Denote 


^  A  d^f{y  = 


-1) 


dy^ 


w{Xi). 


(3.30) 


By  pulling  out  the  complex  time  coefficients  and  denoting  them  as  the  vector  a  as 
above  we  may  construct  a  complex  observation  matrix,  O 


z{xi,y  =  -l,t)  =  Oa  +  Dq{t). 


(3.31) 


Finally,  by  creating  p  by  stacking  the  real  and  imaginary  parts  of  a,  as  well  as  q{t), 
we  may  construct  an  observation  equation  in  state-variable  form 


ZK(Xi>J'  =  -1.0 

Or  —0/  D 

_  zi{xuy  =  -ht) 

0/  Or  0 

(3.32) 


Since  we  may  measure  only  real  output  (shear),  we  are  left  with  only  the  top  half  of 
the  observation  matrix. 


ZRixi,  y  =  -1,  t)  =  [  Or  -Oi  D  ]  p.  (3.33) 

In  order  to  describe  the  system  in  traditional  control  terms,  define 

C^[Or  -Oj  D].  (3.34) 


Then, 


ZR{Xi,y  =  -l,t)  =  Cp 


(3.35) 


3.3,4.  Initial  conditions 

We  have  not  accounted  for  the  initial  value  (2.11)  in  our  boundary  value  problem. 
We  may  account  for  the  initial  condition  by  assuming  it  can  be  written  as  a  series 
expansion  in  terms  of  Pn{x)  and  Then 

<f){x,  y,t  =  0)  =  ip(x,  y,  t  =  0)  -  s{t  =  0)/(y)w(x)  =  g(x,  y)  -  s(t  =  0)/(y)w(x) 

N  M 

=  E  'EbMx)r„{y) 


(3.36) 
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where  b„m  are  assumed  known.  Then 

a{t  =  0)  =  b.  (3.37) 

Since  a  and  b  are  both  complex,  stack  the  real  and  imaginary  components  to  arrive 
at  an  initial  condition  on  p : 


Kt  =  0)  = 


h 

q{t  =  0) 


(3.38) 


We  have  now  defined  the  system  completely  in  the  state-variable  form  of  (3.7),  (3.9) 
where  A  and  B  are  given  in  (3.28),  C  is  given  in  (3.35),  and  the  initial  condition  is 
(3.38).  Although  we  have  shown  the  state-space  formulation  using  a  single  input  and  a 
single  output,  we  may  formulate  the  multiple  input/multiple  output  case  similarly.  In 
that  case,  the  multiple  input  functions  are  represented  as  a  sum  of  separable  functions 
in  (2.12)  and  multiple  outputs  are  represented  as  a  vector  whose  components  are  of 
the  form  (3.35).  It  should  be  noted  that  in  the  multiple  input  case,  it  is  important  for 
physical  implementation  that  for  any  finite  stretch  in  the  x-direction,  only  one  of  the 
multiple  input  functions  is  non-zero. 


4.  Infinite-dimensional  nature  of  the  channel  problem  and  effect  on 
actuator  design 

In  a  system  described  by  partial-differential  equations,  such  as  the  channel  flow 
problem,  no  finite-dimensional  model  will  capture  all  the  dynamics  of  the  system. 
Such  a  system  is  called  infinite-dimensional.  By  examining  the  Galerkin  approximate 
solution  (3.12),  we  see  that  the  approximate  solution  can  converge  only  &sN,M  -*  oo. 
This  would  result  in  an  infinite  number  of  ordinary  differential  equations.  Each  n 
value  included  in  the  approximate  solution  is  referred  to  as  adding  an  additional 
wavenumber  in  the  dynamic  model.  Even  for  just  one  wavenumber,  n  =  N  (say), 
we  see  the  assumed  Galerkin  solution  requires  an  infinite  number  of  terms  in  the 
y-direction, 

M-^00 

Ux,y,  aNn,{t)PN{x)rM  (4.1) 

m=0 

Interestingly,  by  examining  the  structure  of  the  A  matrix  in  the  resulting  state-space 
model,  we  see  that  the  dynamics  associated  with  each  individual  wavenumber  are 
decoupled  from  the  others.  This  is  characterized  by  a  block  diagonal  form  of  the  A 
matrix,  (4.2).  As  a  result,  we  may  separate  the  problem  of  determining  system  poles 
into  a  set  of  smaller  problems  that  include  only  one  wavenumber  at  a  time: 


[4(n=l)]  0  •  •  •  00 

0  [4(„=2)]  0  0 

0  0  0 

00  0  0  [4(„=«,)]  . 


(4.2) 


Indeed,  the  eigenvalues  (poles)  of  the  entire  A  matrix  are  simply  the  eigenvalues 
(poles)  of  each  smaller  A  matrix  block.  Offsetting  this  computational  advantage, 
however,  is  the  fact  that  the  locations  of  the  zeros  are  not  decoupled  by  wavenumber. 
Therefore,  even  though  the  poles  of  the  system  can  be  calculated  separately  using 
only  one  wavenumber  in  a  particular  computation,  the  zeros  of  the  system  force  all 
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the  wavenumbers  in  the  model  to  be  considered  together.  This  can  easily  be  seen 
from  an  examination  of  the  transfer  functions  associated  with  just  two  wavenumbers. 
Consider  two  transfer  functions  obtained  by  using  two  distinct  wavenumbers  and  the 
same  input,  u{t),  and  output,  z{t): 


N(„=\){s) 

'D(„=1)(S)’ 


„  _  ^(n=2)is) 

Di„=2){sy 


(4.3) 


The  roots  of  N(s)  are  the  zeros  and  the  roots  of  D(s)  are  the  poles.  In  terms  of  input 
and  output, 

Z  =  Z  =  H(n=2)U.  (4.4) 

Considered  together,  however. 


z 


[H(n=l)  +  H(„=2)]  M  = 


D(n=2){s)N(„=\){s)  +  D(n=l){s)N{„=2){s) 

D(„=i){s)D(„=2)is) 


(4.5) 


The  poles  of  the  new,  combined  transfer  function  are  clearly  those  of  each  model 
separately.  The  zeros,  however,  are  the  roots  of  a  completely  new  polynomial 
D{„=2)is)N(n=i){s)  +  D(„^i)(s)N(„=2){s),  which,  in  general,  has  Uttle  to  do  with  either  of 
the  original  numerators. 

We  have  seen  that  an  infinite  expansion  is  needed  to  fully  model  the  channel  system. 
Since  in  practice  an  infinite  expansion  cannot  be  used  to  create  a  state-space  model,  we 
will  always  have  some  part  of  the  system  dynamics  that  is  unmodelled.  Furthermore, 
we  have  seen  that  the  A  (dynamic)  matrix  may  be  decoupled  by  wavenumber. 

Consider  a  partition  of  the  state-space  model  as 


r 

dt 

L  dt  J 


Am  0 
0  Au 


Xu 


+ 


Bu 


«(<). 


(4.6) 


z=[Cm  Cu]x,  (4.7) 

where  the  subscripts  m  and  u  represent  the  modelled  and  unmodelled  parts  of  the 
system.  In  the  channel  flow  problem,  the  unmodelled  part  is  meant  to  denote  only 
the  dynamics  of  wavenumbers  left  out  of  the  finite-dimensional  model.  Note  that 
both  Au  and  A„  are  of  infinite  dimension;  Au  because  of  the  infinite  number  of 
wavenumbers  left  out  of  the  reduced-order  model  and  Am  because  of  the  infinite 
number  of  expansion  functions  needed  in  y  for  each  of  the  finite  number  of  modelled 
wavenumbers.  One  way  to  avoid  considering  unmodelled  wavenumber  dynamics  is 
to  ensure  that  the  control  input,  u{t),  has  no  effect  whatsoever  on  the  unmodelled 
wavenumber  dynamics.  In  terms  of  the  state-space  model  (4.6),  this  is  equivalent  to 
rendering  =  0.  This  is  known  as  making  the  unmodelled  wavenumber  dynamics 
uncontrollable.  By  examining  (3.27),  we  see  that  the  B  matrix  is  formed  from  terms 
of  in  (3.16).  If  w(x)  in  (3.16)  is  such  that  =  0,  then  those  components  of  the 
B  matrix  are  zero.  Equivalently,  we  must  ensure  that  the  projection  of  w(x)  onto  the 
unmodelled  wavenumbers  is  zero.  Owing  to  the  orthogonality  of  Fourier  components, 
we  select  w(x)  to  be  made  up  of  modelled  Fourier  components  only. 


5.  Single-wavenumber  channel  model 

Consider  the  channel  model  shown  in  figure  1.  The  total  non-dimensional  length  of 
the  channel,  L,  is  4n.  In  this  case,  the  fundamental  wavenumber,  oo  =  0.5.  Only  one 
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•-dw(x)/dbc  =  -cos(x)  Shear  sensor  at  x  -  +n/2 

Figure  L  System  model  for  Poiseuille  flow.  Input  is  applied  along  the  bottom  plate  as 
w(x)  =  sin(x)  and  shear  is  measured  at  x  =  ;c/2  with  Re  =  10000  and  a  =  1.0. 


wavenumber  is  included  in  the  model,  corresponding  to  a  =  /oq  =  1.00.  The  Reynolds 
number  is  chosen  as  Re  =  10000.  Input  is  distributed  along  the  bottom  plate 
with  a  sinusoidal  weighting  function  in  order  to  render  the  unmodelled  wavenumber 
dynamics  uncontrollable  as  described  in  §4.  In  terms  of  the  input  function,  w(x),  see 
(2.12) 

w(x)  =  sin(x).  (5.1) 

Note  that  the  physical  blowing/suction,  v{x,  y,  t),  takes  the  form 

o(x,y,  t)  =  =  -1)  =  -9(0cos(x)/(y  =  -1)  (5.2) 

see  (2.21).  This  type  of  input  may  be  achieved  in  practice  by  a  large  number  of 
independently  controllable  actuators  that  are  distributed  along  the  lower  channel 
wall.  The  f(y)  function  is  chosen  as  in  (2.20).  The  sensor  location  is  x  =  +\n.  In 
order  to  visualize  the  control-theoretic  model,  the  system  A,  B,  and  C  matrices  are 
transformed  to  transfer  function  form,  H(s).  Figure  2  shows  the  locations  of  the  poles 
and  zeros  in  the  s-plane  for  the  single-wavenumber  model. 


5.1.  Relation  of  poles  to  eigenvalues  of  the  Orr-Sommerfeld  equation 

Poles  of  the  transfer  function  or,  equivalently,  eigenvalues  of  the  A  matrix  are  closely 
related  to  the  eigenvalues  of  the  Orr-Sommerfeld  equation.  Assume  a  solution  of 
(2.6)  of  the  form 

V(x.y,t)^j8(y)e‘«“e-‘“'.  (5.3) 

By  substituting  into  (2.6),  we  obtain  the  familiar  Orr-Sommerfeld  equation  in  the 
normal-mode  form. 


dy^ 


1 


iaRe  [  dy^ 


+  «'‘^(y) 


(5.4) 


Here,  the  Reynolds  number.  Re,  is  known;  the  wavenumber,  a,  is  assumed  real  and 
known;  the  complex  wave  speed,  c,  is  the  eigenvalue  of  the  problem;  and  the  function 
P{y)  is  the  eigenvector  of  the  problem.  Stability  of  a  flow  for  a  given  value  of  Re 
and  a  is  determined  by  the  imaginary  part  of  c.  If  the  imaginary  part  is  positive,  the 
solution  (5.3)  becomes  an  unbounded  exponential  and  flow  is  unstable.  Poles  6f  the 
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Rgure  2.  Pole  (x)  and  zero  (o)  configuration  for  the  system  model  of  figure  1.  Note  that  all  poles 
and  zeros  appear  in  complex-conjugate  pairs.  Each  pair  of  complex-conjugate  poles  represents  one 
mode  of  the  system.  One  pair  of  poles  is  to  the  right  of  the  imaginary  axis.  This  indicates  the  system 
has  one  unstable  mode.  Note  that  the  single  pole  at  the  origin  represents  a  built-in  integrator  due 
to  «(t)  =  dq(t)/dt.  Channel  model:  Re  =  10000,  shear  sensor  at  7t/2,w(x)  =  sin(x),  L  =  47t,  a  =  1.0. 


transfer  function  in  the  Laplace  domain  at  s  =  p,-  correspond  to  solutions  in  the  time 
domain  of  e^'*.  As  a  result,  we  observe  that  poles  of  the  transfer  function  are  related 
to  eigenvalues  of  the  Orr-Sommerfeld  equation  as 

Pi  =  — iac.  (5.5) 

In  order  to  validate  our  linear  code,  we  compare  our  poles  to  eigenvalues  produced  in 
Orszag  (1971).  Orszag  (1971)  obtained  Orr-Sommerfeld  eigenvalues  for  the  channel 
problem  with  a  =  1.00  and  Re  =  10000.  He  reported  only  one  slightly  unstable 
eigenvalue  at  s  =  0.00373967  -1-  iO.23752649  (c  =  0.23752649  -  iO.00373967).  The 
eigenvalue  is  seen  as  unstable  by  its  positive  real  part.  We  obtained  identical  results 
including  the  one  unstable  inode  at 

s  =  0.00373967  T  iO.23752649. 

All  other  stable  modes  obtained  in  the  present  study  are  identical  to  those  reported 
in  Orszag  (1971).  The  goal  of  our  control  system  will  be  to  move  these  unstable  poles 
into  the  stable  half  of  the  s-plane  or,  equivalently,  make  sure  the  controlled-system 
poles  all  have  real  parts  less  than  zero. 

5.2.  Verification  of  model  zeros 

Verification  of  the  system  zeros  is  more  difficult  due  to  the  fact  that  no  published 
results  exist  to  our  knowledge.  We  use  the  channel  code  simulation  used  in  Kim, 
Moin  &  Moser  (1987)  to  verify  our  zeros.  This  code  is  a  spectral  channel  flow  code 
which  uses  periodic  boundary  conditions.  Consider  the  transfer  function  model  (3.2) 
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figure  3.  Output  from  the  Navier-Stokes  simulation  with  input  applied  as  «(t)  = 

Although  an  unbounded  input  is  applied,  near  zero  output  is  observed  at  sensor  location  n/2  (solid 
curve).  Sensor  location  0  (dashed  curve)  and  sensor  location  n  (dash-dot  curve)  do  not  show  zero 
output.  Channel  model:  Re  =  10000,  shear  sensor  at  7t/2,w(x)  =  sin(x),  L  =  43t,  a  =  1.0. 


with  J  =  1  and  /  =  1 : 


In  the  Laplace  domain, 


H{s)  = 


Z(s)  _  (s-Ci) 
Uis)  {s-piY 


sZ{s)-pxZis)  =  sU{s)-t:iUis). 


(5.6) 

(5.7) 


Assuming  zero  initial  conditions  and  transforming  into  the  time  domain,  (5.7)  becomes 

^-piz(t)  =  ^-Ciu(t)  (5.8) 

where  z(f)  is  the  output  function  and  u(r)  is  the  input  function.  If  the  input  is  taken 
as 

u{t)  =  e^'*  (5.9) 

then  the  right-hand  side  of  (5.8)  becomes  zero  and  the  system  behaves  as  if  zero 
input  has  been  applied.  As  a  result,  z{t)  remains  zero  for  all  time.  Therefore,  we  naay 
verify  zeros  of  the  system  by  applying  non-zero  input  in  special  ways  and  observing 
zero  output.  Note  from  figure  2  that  for  sensor  location  x  =  +^n,  one  zero  exists 
in  the  right  half  plane  at  s  =  -FO.317557  -|-  Oi.  This  represents  an  ideal  zero  to  check 
as  this  corresponds  to  an  unbounded,  unstable  input.  Figure  3  shows  output  from 
the  Navier-Stokes  simulation  with  input  u{t)  =  dq{t)/ dt  —  e"*"®-^*’**’*.  Note  that  even 
though  an  unbounded,  exponentially  growing  input  is  applied  to  the  system,  near 
zero  output  is  observed  at  the  sensor,  thus  verifying  that  is  indeed  a  zero  of  the 
system.  This  is  contrasted  with  measurements  at  different  sensor  locations  that  grow 
rapidly.  With  a  sensor  at  a  different  location,  the  zeros  of  the  system  change  position 
so  zero  output  is  no  longer  expected  for  this  particular  input.  Absolute  zero  is  not 
observed  at  sensor  location  x  =  \n  due  to  slight  numerical  inaccuracies  in  the  model. 


6.  Feedback  control  and  stabilization 
At  this  point,  with  the  construction  of  a  valid  state-variable  model,  any  number  of 
control  schemes  may  be  employed  to  stabilize  the  system.  A  general  control  system  is 
shown  in  figure  4.  The  output  of  the  system  is  fed  into  a  controller.  The  output  of  the 
controller  is  then  used  to  create  an  input  to  the  system.  The  design  of  a  controller  that 
achieves  certain  system  characteristics  is  the  goal  of  control  system  design.  Several 
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modem  control  techniques  may  be  applied  that  require  a  state-variable  model.  In  this 
paper,  the  primary  goal  will  be  system  stability.  For  this  purpose,  a  simple  constant 
gain  feedback  with  integral  compensator  will  be  shown  to  be  suflScient. 


6.1.  Constant  gain  feedback  with  integral  compensator 

Figure  5  shows  an  integral  compensator  feedback  control  scheme.  K  is  referred  to  as 
the  gain  of  the  feedback.  The  output  of  the  system,  in  this  case  shear,  is  multiplied 
by  a  feedback  gain,  integrated  in  time,  and  then  fed  back  as  blowing  and  suction  at 
the  input.  Note  that  the  integrator  is  required  because  of  the  stracture  of  the  input 
given  in  (3.27),  i.e.  the  input  is  taken  as  the  derivative  of  suction/blowing.  The  signal, 
r(t),  is  called  the  reference  signal.  It  is  used  to  define  the  desired  output.  In  our  case, 
we  would  like  the  shear  output,  z(x„y  =  — l,t),  to  be  zero.  Therefore,  we  set  r(t)  to 
zero.  Then 

u(t)  =  =  -Kz{Xi,  y  =  -1, 0-  (6-1) 

Therefore,  assuming  q{0)  =  0, 

q{t)  =  -K  [  z{Xi,y  ^ -U'c)dz  (62) 

Jo 

and  by  observing  (2.21),  we  may  describe  physical  blowing  and  suction  at  the  bound¬ 
ary: 

?(x,y  =  -1,0  =  =  -l)  =  K cos(x)  jT  z(xi,y  =  -l,T)dT  (6.3) 


As  defined  earlier, 

i?[u(0]  UisY 


(6.4) 


Therefore, 


Z(s)  =  His)U{s).  (6.5) 

In  the  absence  of  feedback,  one  mode  is  unstable.  Also,  one  pole  exists  at  the  origin 
for  the  integrator.  Consider  a  new  transfer  function  from  the  reference  input,  r(f),  to 
the  output,  z(t),  in  the  presence  of  feedback: 


u{t)  =  r{t)  -  Kz{t). 


(6.6) 


By  taking  the  Laplace  transform  of  (6.6), 


U{s)  =  R{s)-KZ{s). 


(6.7) 


Then 


Finally, 


Z(s)  =  His)[Ris)-KZ{s] 

Zjs)  Hjs) 

R{s)  l-KH{sy 

The  new  poles  of  the  feedback  system  are  defined  by 

l-KH(s)  =  0. 


(6.8) 

(6.9) 

(6.10) 


As  K  gets  larger  and  larger,  it  is  clear  that  the  poles  of  the  new  system  tend  toward 
the  zeros  of  H{s).  In  this  way,  modes  of  the  system  can  sometimes  be  changed  to 
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Figure  4.  Feedback  control.  The  output  of  the  system  is  fed  into  a  controller  and  then  to  the  input. 


Feedback  gain 


Rgure  5.  Feedback  control  for  channel  system.  The  output  of  the  system  is  multiplied  by  a 
feedback  gain,  integrated  in  time,  and  then  fed  back  at  the  input  The  reference  signal,  r(t), 
equals  0. 

fonn  a  stable  system  (Franklin  et  al.  1988).  Unstable  modes  that  appear  as  poles 
on  the  right  hand  side  of  the  complex  s-plane  and  ‘marginally  stable’  poles  on  the 
Im(s)-axis  are  drawn  to  the  left-hand  side  by  applying  feedback. 

6.2.  Sensor  placement 

We  have  seen  that,  by  applying  feedback,  poles  of  the  system  will  eventually  be  drawn 
to  zeros  of  the  system.  In  the  channel  flow  system  of  flgure  2,  any  feedback  will 
cause  either  the  pole  at  the  origin  (integrator  pole)  or  the  unstable  mode  poles  to  be 
drawn  to  the  zero  on  the  real  axis  in  the  right-hand  plane,  thus  making  the  system 
unstable  to  a  greater  degree.  Therefore,  finding  a  transfer  function  that  has  all  zeros 
in  the  left-hand  plane  becomes  an  important  objective.  The  poles  of  the  system  are 
independent  of  sensing  or  actuation.  However,  the  zeros  of  the  system  are  dependent 
on  both  the  type  and  location  of  sensing  and  the  type  and  location  of  actuation. 
Figure  6  shows  the  pole/zero  configuration  for  the  channel  model  with  the  shear 
sensor  at  three  different  locations.  Only  the  top  half  of  the  s-plane  is  shown  since 
the  bottom  half  is  a  mirror  image  projected  across  the  real  axis  as  shown  in  figure  2. 
The  poles  of  all  the  models  are  in  the  same  location  as  expected.  However,  the  zeros 
are  different  in  all  three  cases.  In  figures  6(a)  and  6(b),  we  observe  a  lone  zero  in 
the  right-hand  s-plane.  However,  when  the  sensor  is  placed  at  x  =  -f-Ti,  figure  6(c) 
shows  all  zeros  in  the  left-hand  plane.  In  fact,  there  is  a  region  around  x  =  n  that 
results  in  all  zeros  in  the  left-hand  plane.  By  placing  a  shear  sensor  at  x  =  rc,  simple 
feedback  with  integral  compensation  will  allow  stabilization  with  the  proper  value 
of  gain,  K.  In  the  case  of  sensor  locations  that  result  in  right-hand-plane,  so  called 
‘non-minimum  phase’,  zeros  stabilization  is  still  possible.  However,  more  complex 
controllers  (Bryson  &  Ho  1975;  Ogata  1990)  must  be  designed  that  are  beyond  the 
scope  of  this  paper. 

6.3.  Root  locus  analysis  and  numerical  simulation 
One  way  to  visualize  how  system  poles  will  change  as  the  feedback  gain,  K,  changes 
is  to  construct  a  root  locus  plot.  This  is  a  plot  of  all  poles  of  a  system  as  the  feedback 


Res 

figure  6.  Pole  (x)/zero  (o)  configuration.  Channel  model:  Re  =  10000,  w(x)  =  sin(x),  L  =  4jt, 
a  =  1.0.  Only  the  top  half  of  the  s-plane  is  shown.  Shear  sensor  at  (a)  jc/4  (b)  3jt/4,  (c)  n. 

gain  varies  from  K  =  0  to  K  =  co.  Figure  7  shows  such  a  plot  for  the  system  shown  in 
figure  6(c).  When  K  reaches  0.1,  all  feedback  system  poles  (closed-loop-poles)  Ue  on 
the  left-hand  plane  and  no  instabilities  exist  in  the  new  feedback  system.  Numerical 
results  obtained  from  the  Navier-Stokes  simulation  for  the  new  feedback  controlled 
system  are  shown  in  figure  8.  The  computation  is  carried  out  without  feedback  until 
t  =  50,  after  which  the  feedback  is  turned  on.  We  see  that  the  growing  instability  is 
quickly  suppressed.  At  the  instant  the  controller  is  turned  on,  the  simulation  shows 
a  high  transient  response  due  to  the  non-continuous  nature  of  the  input  at  that  time 
instant.  In  terms  of  the  state  space  defined  in  (3.28),  (3.35), 


(6.11) 

zixi,y  =  -l,t)  =  Cp, 

(6.12) 

u{t)  =  -Kzixi,y  =  -l,t). 

(6.13) 

Then,  in  a  closed  loop. 

^=Ap-  B{Kz{xu  y  =  -l,t)) 

(6.14) 

^Ap-  BKCp  =  (4  -  BKC)p 

(6.15) 
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Figure  7.  Root  locus  plot  for  channel  system  with  input  as  in  the  channel  model  (figure  1)  and 
shear  output  at  x  =  te.  The  poles  start  at  the  open-loop  poles  shown  with  x’s.  As  K  increases,  they 
start  to  move  toward  the  location  of  the  system  zeros,  shown  as  o’s.  The  pole  at  position  s  =(0,0) 
moves  directly  to  the  left.  The  unstable  pole  moves  quickly  to  a  position  just  to  the  left  of  the 
imaginary  axis.  Near  0.7  <  y  <  1.0,  -’0.6  <  x  <  —0.4,  we  see  a  pole  moving  towards  a  zero  that 
is  out  of  the  range  of  this  figure.  Channel  model:  Re  =  10000,  shear  sensor  at  te,  w(x)  =  sin(x), 
L  =  47e,  a  =  1.0.  Only  the  top  half  of  the  s-plane  is  shown. 


figure  8.  Navier-Stokes  simulation  of  feedback  control.  Shear  is  multiplied  by  a  feedback  gain, 
integrated,  and  fed  back  into  the  input.  Channel  model:  Re  =  10000,  feedback  shear  sensor  at  te, 
w(x)  =  sin(x),  L  =  47e,  a  =  1.0.  Curve  1  (solid)  shows  shear  at  location  te;  curve  2  (dot-dash)  shows 
shear  at  location  7e/2. 

The  solution  to  this  equation  is 

(6.16) 

This  solution  will  approach  zero  as  t  ^  oo  since  the  eigenvalues  of  4  —  BKC  (closed- 
loop  poles)  are  all  stable. 

6.4.  Robustness  in  the  presence  of  Reynolds  number  uncertainty 

A  major  advantage  of  feedback  control  systems  is  their  robustness  to  system  un¬ 
certainty.  From  a  practical  point  of  view,  Reynolds  numbers  may  not  be  known 
exactly  or  may  change  frequently  as  in  the  flight  of  an  airplane,  for  example. 
Figure  9  shows  the  open-loop  pole/zero  configurations  for  systems  with  varying 
Reynolds  numbers.  Note  that  these  systems  start  with  all  zeros  in  the  left-hand 
s-plane.  A  root  locus  analysis  shows  that  a  feedback  system  with  X  =  0.1  stabilizes 
both  systems.  Indeed,  a  feedback  gain  of  K  =  0.1  stabilizes  systems  for  a  wide 
range  of  Reynolds  numbers  from  Re  =  1000  to  Re  =  40000.  Figure  10  shows 
the  least-stable  pole  in  both  the  controlled  and  uncontrolled  systems  for  several 
Reynolds  numbers.  Recall  that  a  pole  with  real  part  less  than  zero  is  stable.  We 
see  near  Re  —  5772,  an  unstable  pole  appears  in  the  open-loop  (uncontrolled) 
system.  Unstable  eigenvalues  continue  to  exist  in  the  open-loop  system  until 
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Figure  9.  Pole  (x)/zero  (o)  configuration.  Channel  model:  shear  sensor  at  n,  w(x)  ==  sin(x), 
L  =  47r,  a  =  LO.  Only  the  top  half  of  the  s-plane  is  shown,  (a)  Re  —  7500,  (b)  Re  =  12  500. 


Rgure  10.  Real  part  of  least-stable  pole  for  open-loop  (solid  line)  and  closed-loop  (dashed  line) 
channel  system  with  1000  <  Re  ^  40000.  Shear  sensor  at  tc,  w(x)  =  sin(x),  L  —  ct  —  1.0.  Data 
points  are  shown  as  x’s. 


30000  <  Re  ^  35000.  An  identical  feedback  controller,  with  gain  R  =  0.1, 
however,  stabilizes  the  system  for  Reynolds  numbers  in  the  range  1000  ^  Re  ^ 
40000.  Clearly,  the  feedback  controller  is  extremely  robust  to  changes  in  Reynolds 
number. 
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7.  Unobservable  transient  response  as  a  path  of  bypass  transition 

Recently,  some  authors  (Trefethen  et  al  1993;  Butler  &  Farrell  1992;  Farrell  1988; 
Henningson  1994)  have  suggested  a  possible  path  of  bypass  transition  in  Poiseuille 
flow  that  is  caused  by  high  transient  response  due  to  the  non-self-adjointness  of  the 
evolution  matrix.  This  high-transient  behaviour  is  due  to  the  non-orthogonality  of  the 
eigenmodes  and  not  caused  by  any  single  mode,  but  rather  a  combination  of  many 
modes.  In  this  section,  we  suggest  a  possible  different  form  of  bypass  transition  in  a 
controlled  system  due  to  a  single  mode.  In  the  presence  of  any  type  of  control,  it  is 
possible  that  energy  from  an  input  is  fed  into  a  mode  such  that  the  mode  is  reinforced 
during  the  transient  period  before  it  eventually  dies  out.  During  the  transient  period, 
however,  the  mode  may  grow  to  an  amplitude  large  enough  to  trigger  nonlinear 
effects  that  induce  transition  to  turbulence.  This  may  be  a  possible  path  of  bypass 
transition  in  a  controlled  system. 

Single  modes  may  be  characterized  in  linear  systems  in  terms  of  modal  controlla¬ 
bility  and  modal  observability.  Modal  controllability  implies  that  a  particular  mode 
may  be  affected  by  the  actuators  chosen  (blowing/suction  in  our  case).  Modal  ob¬ 
servability  implies  that  a  particular  mode  may  be  measured  by  the  sensors  chosen 
(shear  in  our  case).  A  particular  mode  may  be  non-observable,  non-controllable,  or 
both.  This  will  be  seen  in  §§7.1  and  7.2.  If  high-transient  modes  are  controllable  and 
observable,  control  theory  may  be  used  to  suppress  them.  However,  if  high-transient 
modes  are  unobservable,  the  high  transients  will  never  be  seen  and  feedback  control 
cannot  be  used  to  suppress  them. 

7.1.  Modal  canonical  form 

The  state-space  formulation  provides  an  excellent  framework  for  assessing  the  rein¬ 
forcement  of  each  mode  individually.  Consider  an  n-dimensional  state  space  with 
scalar  input  and  output 

^  =  Ax  +  Bu,  (7.1) 

at 

z  =  Cx.  (7.2) 

We  may  perform  a  similarity  transformation  on  the  system  to  produce  a  new  state- 
space  representation  with  the  same  input-output  relationship,  but  with  a  different 
interpretation  of  the  state  variables.  Let 

P  =  [  vi  V2  ...  v„  ]  (7.3) 

where  vj  is  the  jth  eigenvector  of  the  A  matrix  and  n  is  the  dimension  of  the  A  matrix. 
A  new  representation  is  constructed  as  (Grace  et  al.  1992;  Kailath  1980) 

—  =  p-^APx  -f-  p-^Bu,  (7.4) 

df 

z  =  CPx:  (7.5) 

where  x  =  P~^x.  The  new  representation  is  written  as 

^  =  A5c  +  Bu,  (7.6) 

dt 

z  =  Cx  (7.7) 

where  the  real  eigenvalues  of  the  original  A  matrix  appear  on  the  diagonal  of  A  and 
the  complex  eigenvalues  appear  in  a  2  x  2  block  on  the  diagonal  of  A,  For  an  A 
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matrix  with  eigenvalues  (aijC  +j(o,(x.2),  the  A  matrix  is 


■  ai  0  0  0  ' 

0  a  CO  0 

0  —ft)  ff  0 

_  0  0  0  a2  . 


(7.8) 


In  this  form,  each  state-variable  pair  represents  a  mode  of  the  system.  Furthermore, 
the  modes  are  block  decoupled  so  that  each  mode  is  represented  by  the  2  x  2  system 


(7.9) 

(7.10) 


where  co,  a,  bi,  h,  ci,  and  C2  are  all  scalars.  This  is  known  as  the  modal  canonical 
state-space  form.  All  modes  can  be  monitored  directly  for  high-transient  reinforcement 
in  this  form. 


7.2.  Modal  observability  and  controllability 

It  is  easy  to  see  that  if  hi  and  62  are  zero,  no  input  can  affect  the  mode  and  the 
mode  is  uncontrollable.  Similarly,  if  ci  and  C2  are  zero,  then  the  motion  of  xi  and 
X2  cannot  be  measured  at  the  output,  zi,  and  the  mode  is  unobservable.  In  terms 
of  stabilization,  an  uncontrollable  mode  is  not  problematic  as  long  as  it  is  stable. 
However,  if  an  unstable  mode  is  uncontrollable,  nothing  can  be  done  to  stabilize  it. 
Observability  may  be  a  problem  even  if  modes  are  stable.  If  an  unobservable  mode 
is  highly  amplified  by  control  energy,  no  attempt  could  be  made  to  suppress  it  since 
it  would  not  be  observed  at  the  output.  In  terms  of  poles  and  zeros,  uncontrollable 
or  unobservable  modes  both  show  up  as  pole/zero  cancellations  in  the  complex  s- 
plane.  As  can  be  seen  from  figure  2,  many  poles  and  zeros  for  the  channel  system 
lie  on  top  of  each  other,  indicating  that  certain  modes  in  the  system  are  either 
uncontrollable,  unobservable,  or  both.  Although  a  mode  may  not  be  physically 
observable  at  the  output,  the  modal  canonical  formulation  allows  us  to  numerically 
observe  the  state  evolution  of  each  mode  directly.  In  this  way,  we  may  assess  the  risk 
of  highly  amplified  modes  triggering  bypass  transition.  This  is  done  for  the  controlled 
system  simulated  in  figure  8.  The  uncontrolled  linear  system  is  simulated  with  an 
initial  condition  for  200  time  steps  at  which  time  a  feedback  controller  with  gain 
K  =  0.1  is  turned  on.  Figure  11  shows  the  state  evolution  of  the  modes  with  poles 
at  s  =  —0.1474  ±  0.8514i  and  s  =  -0.3252  ±  0.6361i  as  well  as  the  state  evolution  of 
the  one  unstable  mode.f  In  addition,  figure  11(a)  shows  the  shear  measurement  at 
the  output.  We  see  that  the  unstable  mode  dies  out  quickly  as  soon  as  the  feedback 
is  turned  on.  In  addition,  the  mode  at  s  =  -0.1474  ±  0.85 14i  gains  almost  no  energy 
after  the  feedback  is  activated.  This  indicates  that  the  mode  is  not  a  bypass  mode. 
The  mode  at  s  =  -0.3252  ±0.6361i  is  seen  to  gain  a  lot  of  energy  after  the  controller 
is  started.  Indeed,  the  amplitude  of  the  transient  response  is  nearly  twice  that  of 
the  unstable  mode.  This  represents  a  possible  bypass  mode  since  such  an  amplitude 

t  Note  that  amplitudes  shown  as  a  result  of  linear  system  simulations  should  only  be  used  for 
comparison  with  each  other  since  the  linear  model  has  been  scaled  to  unity  open-loop  feedforward 
gain,  i.e.  H(s)  =  icZ(s)/l/(s)  where  »c  =  1. 
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Figure  11.  Linear  model  state  evolution  of  system  modes  with  control  applied  at  t  —  200.  {a) 
measured  shear,  (h)  the  state  evolution  of  the  unstable  mode,  (c)  the  state  evolution  of  the  mode 
at  s  =  -0.1474  +  0.8514i,  (d)  the  state  evolution  of  the  mode  at  s  =  -0.3252  +  0.6361i.  Channel 
model:  Re  —  10000,  shear  sensor  at  n,  w(x)  =  sin(x),  L  =  47r,  a  =  1.0. 


may  force  the  system  into  the  nonlinear  region  where  transition  to  turbulence  may  be 
triggered.  Furthermore,  the  high  energy  nature  of  the  mode  cannot  be  observed  at  the 
shear  sensor  since  this  mode  is  unobservable.  This  can  be  seen  from  the  plot  of  the 
shear  sensor  output.  Modal  un-observability  is  also  suggested  by  noting  the  relatively 
low  Cl  and  C2  values  of  0.5155,  -0.1131  compared  to  the  ci  and  C2  values  of -162.85, 
234.66  for  the  unstable  mode  which  is  clearly  seen  at  the  output.  Fortunately,  the  Wgh 
energy  nature  of  the  mode  does  not  lead  to  a  bypass  transition  in  this  case.  This  is 
verified  by  the  Navier-Stokes  simulation  in  figure  8,  which  shows  no  nonlinear  effects. 
If  the  instability  were  allowed  to  grow  to  a  higher  amplitude  before  the  controller 
was  applied,  however,  the  highly  amplified  mode  in  the  transient  response  might  have 
triggered  nonlinear  effects. 


8.  Multiple  instability  control 

Although  most  work  has  focused  on  suppression  of  single  instabilities,  more  realistic 
models  should  include  multiple  wavenumbers.  Indeed,  for  a  given  Reynolds  number, 
an  infinite  number  of  wavenumbers  exist.  Each  wavenumber  contributes  its  own  poles 
and  zeros  to  the  control-theoretic  model.  Consider  a  model  with  non-dimensional 
channel  length  207r,  where  input  is  applied  as  w(x)  =  sin(x)  +  sin(0.9x).  In  this  model, 
wavenumbers  of  0.9  and  1.0  are  included.  Both  wavenumbers  lead  to  unstable  modes. 
The  pole/zero  configuration  of  this  new  two-wavenumber  model  with  shear  sensor 
at  n  is  shown  in  figure  12(a).  Note  that  in  comparison  to  the  one-wavenumber 
model,  more  poles  and  zeros  exist.  The  original  one-wavenumber  model  contained 
one  *fork’  structure  of  poles,  while  the  two-wavenumber  model  contains  two  fork 
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Rgure  12.  (o)  Pole  (x)  and  zero  (o)  configuration  of  channel  system  of  length  207t,  including 
wavenumbers  of  0.90  and  1.00.  Re  =  10000,  shear  sensor  at  n,  w(x)  =  sin(x)  +  sin(0.9x).  Only  the 
top  half  of  the  s-plane  is  shown,  (b)  Same  as  (a)  but  showing  closed-loop  poles  after  feedback  with 
gain  K  =  0.1. 


figure  13.  Navier-Stokes  simulation  of  channel  system  of  length  207t,  including  wavenumbers  of 
0.90  and  1.00.  Re  =  10000,  feedback  shear  sensor  at  n,  w(x)  =  sin(x)  -b  sin(0.9x),  feedback  gam 
K  =0.1.  Curve  1  (solid)  shows  shear  at  location  it;  curve  2  (dash-dot)  shows  shear  at  location  2n. 


structures.  It  can  be  visualized  that  in  a  model  with  several  wavenumbers,  several 
‘forks’  will  stack  on  top  of  each  other  in  the  s-plane.  Near  s  =  iO.2  are  two 
unstable  poles  in  the  right-hand  plane  that  represent  the  two  unstable  modes  in  the 
system.  As  seen  before,  all  zeros  he  in  the  left  hand  s-plane.  Figure  12(b)  shows 
the  closed-loop  poles  after  feedback  with  gain  K  =  0.1.  Results  from  the  Navier- 
Stokes  simulation  are  shown  in  figure  13.  The  computation  is  carried  out  without 
feedback  until  t  =  200.  A  combination  of  two  grovnng  waves  is  seen  at  the  shear 
sensor.  At  t  =  200,  feedback  with  integral  compensation  is  applied  and  the  system  is 
stabilized. 
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9.  Effect  of  linear  controllers  on  two-dimensional  finite-amplitude 
disturbances 

The  critical  Reynolds  number  above  which  linear  instabilities  exist  in  plane 
Poiseuille  flow  is  Re  «  5772.  Below  this  Reynolds  number,  Poiseuille  flow  is  linearly 
stable.  However,  transition  from  laminar  to  turbulent  flow  occurs  in  experiments 
at  much  lower  Reynolds  numbers,  typically  around  Re  =  1000.  Moreover,  even 
for  super-critical  Reynolds  number  flows,  linear  instabilities  predict  extremely  slow 
growth  rates.  In  practice,  transition  occurs  orders  of  magnitude  more  quickly.  There¬ 
fore,  linear  stability  alone  does  not  dictate  transition.  Transition  to  turbulence  is 
generally  accepted  to  be  a  nonlinear  three-dimensional  phenomenon.  In  previous  sec¬ 
tions,  we  constructed  two-dimensional  linear  controllers  based  on  the  two-dimensional 
lineari2Bd  Navier-Stokes  equations.  By  applying  these  controllers  to  plane  Poiseuille 
flow  with  infinitesimal  two-dimensional  disturbances,  we  saw  that  we  may  linearly 
stabilize  the  system  so  that  the  flow  could  no  longer  support  those  disturbances.  If 
we  apply  our  linear  controllers  to  flows  that  contain  finite-amplitude  disturbances,  we 
can  no  longer  use  linear  analysis  to  describe  the  flow  dynamics  as  nonlinear  terms 
become  relevant.  However,  the  application  of  a  controller  based  on  Unear  analysis 
does  change  the  nonlinear  system. 

Many  authors  (Orszag  &  Patera  1983;  Bayly,  Orszag  &  Herbert  1988)  have  explored 
the  effect  of  two-dimensional  finite-ampUtude  disturbances  on  three-dimensional  in¬ 
finitesimal  disturbances  in  plane  Poiseuille  flow.  Above  Re  «  2900,  for  two  dimensions, 
stable  non-attenuating  finite-ampUtude  equiUbria  exist  in  plane  Poiseuille  flow.  Below 
Re  »  2900,  non-decaying,  finite-ampUtude  equiUbria  do  not  exist.  However,  in  flows 
with  1000  <  Re  <  2900,  the  timescale  for  decay  is  so  large  that  the  flow  may  be 
considered  in  ‘quasi-equiUbrium’.  It  has  been  shown  that  in  the  presence  of  such 
two-dimensional  finite-ampUtude  disturbances,  infinitesimal  three-dimensional  distur¬ 
bances  are  highly  unstable  and  may  cause  transition  in  shear  flows.  Figure  14  shows 
the  energy  of  a  single-wavenumber,  two-dimensional  finite-ampUtude  disturbance 
(a  =  1.0)  and  a  single-wavenumber-pair,  three-dimensional  infinitesimal  disturbance 
(a  =  l.O,  P  =  Tl.O)  obtained  through  direct  numerical  simulation  at  Re  =  3000.  The 
maximum  ampUtude  of  the  two-dimensional  finite-ampUtude  disturbance  is  0.117c. 
At  this  Reynolds  number,  wavenumber,  and  initial  energy  level,  we  see  .the  two- 
dimensional  finite-ampUtude  disturbance  decaying  slowly.  The  three-dimensional  dis¬ 
turbance,  on  the  other  hand,  is  seen  to  rapidly  gain  energy.  Orszag  &  Patera  (1983) 
show  that  the  two-dimensional  instability  acts  as  a  mediator  for  transfer  of  energy 
from  the  mean  flow  to  the  three-dimensional  disturbance,  but  does  not  directly  provide 
energy.  The  growth  rate  of  the  three-dimensional  disturbance  is  orders  of  magnitude 
larger  than  that  of  the  Orr— Sommerfeld  instabiUties. 

Orszag  &  Patera  show,  in  the  case  where  Re  ^  1000,  that  the  attenuation  of 
finite-ampUtude  two-dimensional  disturbances  is  large  enough  that  three-dimensional 
disturbances  do  not  become  unstable.  The  fact  that  high  attenuation  of  the  two- 
dimensional  finite-ampUtude  disturbance  prevented  three-dimensional  instability  be¬ 
low  Re  »  1000  suggests  that  if  controllers  can  be  created  that  speed  up  the  attenuation 
of  the  finite-ampUtude  two-dimensional  disturbance  for  flows  with  Reynolds  num¬ 
bers  greater  than  1000,  three-dimensional  instability  may  be  eliminated.  We  have 
seen  in  §6  that  Unear  controllers  did  stabiUze  infinitesimal  two-dimensional  distur¬ 
bances.  Figure  15  shows  the  effects  of  the  linear  controller  of  §6  when  appUed 
to  a  system  with  the  same  finite-ampUtude  two-dimensional  disturbance  and  in¬ 
finitesimal  three-dimensional  disturbance  shown  in  figure  14.  The  linear  controller 
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Figure  14.  Energy  of  a  two-dimensional  finite-amplitude  disturbance  (a  --  l-O)  ^nd  a 
three-dimensional  infinitesimal  disturbance  (a  =  1.0,^  =  +1.0)  for  Re  =  3000  Solid  line  rep¬ 
resents  the  total  energy  of  the  two-dimensional  finite-ampUtude  disturbance  and  the  dotted  line 
represents  that  of  the  three-dimensional  infinitesimal  disturbance. 


FIGURE  15.  Effect  of  linear  controller  when  applied  at  t  =  15  for  Re  -  3000.  The  *^°  .**““ 
show  the  two-dimensional  finite-amplitude  disturbance  without  (solid)  and  with  (dash-dot)  linear 
control.  We  see  at  t  =  15  that  the  controlled  two-dimensional  finite-amplitude  disturbanee  is  ™g"*y 
damped.  The  bottom  two  lines  show  the  infinitesimal  three-dimensional  disturbance  without  (dot) 
and  with  (dashed)  Unear  control.  That  without  controller  gains  energy  qmckly,  while  that  with 
controller  quickly  fails  to  gain  energy. 
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dramatically  increased  the  attenuation  of  the  finite-amplitude  two-dimensional  dis¬ 
turbance.  As  a  result,  the  infinitesimal  three-dimensional  disturbance  was  rendered 
stable.  Clearly,  the  linear  controller  not  only  stabilized  the  two-dimensional  lin¬ 
ear  system,  but  also  had  a  stabilizing  effect  on  the  three-dimensional  nonlinear 
system.  These  results  show  the  promise  of  linear  controllers  even  in  nonlinear  sys¬ 
tems. 


10.  Conclusions 

In  this  paper,  we  have  developed  feedback  controllers  that  Unearly  stabilize 
plane  Poiseuille  flow.  We  used  a  Galerkin  spectral  method  to  generate  state-space 
models.  Indeed,  any  (convergent)  numerical  method  that  reduces  the  governing  par¬ 
tial  differential  equations  into  a  set  of  ordinary-differential  equations  may  be  used. 
It  is  important  to  note,  however,  that  the  meaning  of  the  state  variables  in  the  state- 
space  model  changes  as  different  numerical  methods  are  employed.  Even  though 
state  variables  may  not  have  any  specific  physical  meaning,  some  numerical  meth¬ 
ods  may  result  in  favourable  divisions  of  system  dynamics.  In  the  plane  Poiseuille 
flow  case,  the  spectral  Galerkin  method  with  the  use  of  Fourier  components  in  the 
x-direction  and  combination-Chebyshev  polynomials  in  the  y-direction  led  to  mod¬ 
elled  dynamics  that  were  de-coupled  by  wavenumber.  This  led  to  a  block  diagonal 
form  of  the  A  matrix.  As  a  result,  we  were  clearly  able  to  describe  modelled  and 
unmodelled  dynamics  in  terms  of  wavenumber  dynamics  included  and  not  included 
in  our  finite-dimensional  model.  This  also  led  us  to  the  concept  of  using  distributed 
control  to  render  certain  wavenumbers  ‘uncontrollable’.  If  other  numerical  methods 
had  been  used,  this  concept  would  not  have  been  so  transparent.  Furthermore, 
plane  Poiseuille  flow  can  be  linearly  stabilized  with  simple  feedback  controllers  if 
sensors  are  placed  at  judicious  locations.  Both  the  position  and  the  type  of  sensing 
and  actuation  change  the  zeros  of  single-input/single-output  models.  In  our  case, 
we  have  seen  that  certain  shear  sensor  locations  lead  to  ‘minimum-phase  systems 
and  some  locations  lead  to  ‘non-minimum  phase’  systems.  Minimum-phase  systems 
are  in  general  easier  to  control  than  non-minimum  phase  systems.  As  a  result,  we 
were  able  to  stabilize  plane  Poiseuille  flow  with  a  constant  gain  feedback,  integral 
compensator  controller.  In  addition,  the  controller  was  extremely  robust  to  a  wide 
range  of  Reynolds  numbers.  Also,  we  have  shown  that  one  danger  of  feedback 
control  is  that  the  linear  transient  response  of  a  controlled  system  can  lead  to  high 
amplitudes  for  a  short  period  of  time.  If  these  amplitudes  are  high  enough,  it  is 
possible  that  they  may  invalidate  the  linear  model  and  enhance  nonlinear  effects. 
Furthermore,  the  high  transients  may  not  be  observable  at  the  output  so  that  feed¬ 
back  control  cannot  be  used  to  suppress  them.  Finally,  we  have  shown  that  linear 
controllers  have  a  strong  stabilizing  effect  on  two-dimensional  finite-amplitude  dis¬ 
turbances.  As  a  result,  three-dimensional  secondary  instabilities  can  be  rendered 
stable. 
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